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PREFACE 


This  document  is  the  final  report  on  the  "Rocket  Plume  Optical  Signature 
Program."  This  work  was  performed  for  General  Research  Corporation, 

1501  Wilson  Boulevard,  Suite  700,  Arlington,  Virgnia  22209,  under  Contract 
040-71-10.  This  contract  covered  the  period  from  20  August  1971  to 
30  September  1972. 

The  effort  described  herein  was  accomplished  by  the  following  study 
personnel:  S.  S.  Cherry,  M.  Thomas,  and  R.  L.  Younkin.  This  report  was 
approved  by  H.  Hurwicz,  Chief  Advance  Technology  Engineer,  Aero/ 
Thermodynamics  and  Nuclear  Effects  Research  and  Development,  Advance 
Systems  and  Technology,  McDonnell  Douglas  Astronautics  Company-West. 


ABSTRACT 


As  part  of  the  Optical  Signatures  Program,  McDonnell  Douglas 
Astronautics  Company-West  has  developed  the  initial  working  model  of  a 
code  to  describe  the  gross  features  of  rocket-plume  radiation  for 
altitudes  above  75  n  mi.  The  main  effort  is  the  construction  of  a 
scheme  for  integration  of  an  arbitrary  function  through  an  arbitrary 
axisymmetric  rocket  plume,  with  any  specified  look  angle,  plume  direction, 
and  vehicle  velocity  direction.  Radiances  are  presented  as  integrated 
values  in  a  specified  spectral  band.  The  equations  used  and  a  printout 
of  the  code  and  of  a  sample  application  are  included. 
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Section  1 

INTRODUCTION  AND  SUMMARY 

A  code  has  been  developed  to  describe  the  gross  features  of  rocket-plume 
radiation  for  altitudes  above  75  mui  (P2 170-FLAME).  The  effort  was 
intended  to  provide  an  initial  working  model,  consistent  '-'ith  available 
technology.  Theoretical  studies  were  not  intended  to  be  conducted  in  any 
of  the  disciplines  required.  The  intent  was,  however,  to  provide  a  suitable 
framework  for  the  description  of  a  plume  with  sufficient  flexibility  to  incor¬ 
porate  the  results  of  future  research. 

This  objective  was  met  by  considering  the  two  most  important  far-field 
radiative  mechanisms,  particle  emission  and  atmospheric  excitation,  in  their 
present  state  of  understanding.  With  these  major  elements  of  plume  radiation 
included  even  in  crude  form,  future  program  modification  to  include  revised 
theories  should  not  require  complete  program  revision.  The  main  effort 
which  has  been  pursued  in  FLAME-code  development  is  the  construction  of 
a  scheme  for  integration  of  an  arbitrary  function  through  an  arbitrary  axi- 
symmetric  rocket  plume,  with  any  specified  look  angle,  plume  direction, 
and  vehicle  velocity  direction.  The  assumption  of  flow  symmetry  about  the 
plume  axis  may  prove  poor  for  late-time  plumes,  but  to  include  a  complete 
three-dimensional  capability  would  have  burdened  the  code's  development  to 
preclude  an  operational  version  at  the  end  of  this  contract.  The  switch  to  3-D 
is  not  difficult  at  a  later  time  and  will  involve  obtaining  flow  correlations  as  a 
function  of  the  azimuthal  angle,  as  well  as  r  and  6,  and  integrating  over  four 
quadrants  instead  of  two. 

Radiances  are  presentee  in  the  current  code  as  integrated  values  in  a  speci¬ 
fied  spectral  band.  Spectral  intensities  are  readily  available  for  the  particles, 
but  only  a  band  theory  is  available  for  the  atmospheric  excited  species,  due 
primarily  to  the  species  cooling  characteristics.  Again,  a  switch  to  spectral 
curves  is  possible  with  ease  at  a  later  time  when  a  more  general  cooling 
theory  evolves. 
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A  schematic  of  the  FLAME  code  is  presented  in  Figure  1-1.  Radiative 
cooling  of  both  the  hot  particles  in  multiphase  plumes  and  the  excited  molec¬ 
ular  radiators  is  included.  The  four  boxes  feeding  into  the  FLAME  code 
represent  the  n  ain  areas  programmed  during  this  contract.  The  use  of  a 
simplified  flow  model  reduces  both  computer  storage  requirements  and  run 
time,  and  allows  a  more  extensive  checkout  of  the  radiative  mechanisms. 
More  sophisticated  flow  models  exist  and  could  be  substituted  at  a  later 
time. 

Input  to  the  code  includes  simple  engine  parameters  (Section  3,  Table  3-1). 
Output  is  a  plot  of  lines  of  constant  radiance  in  the  plume  as  viewed  in  the 
specified  direction.  A  typical  output  from  the  code  is  shown  in  Figure  1-2 
for  a  broadside  look  at  a  solid  propellant  plume  in  the  absence  of  atmospheric 
excitation.  The  radiance  is  plotted  in  a  plane  perpendicular  to  the  direction 
of  view.  The  plume  of  Figure  1-2  is  shown  again  in  Figure  1-3  where  the 
observer  is  at  a  45-degree  angle  to  the  plume  axis.  Radiance  from  the  same 
engine  at  low  altitudes,  including  the  effect  of  atmospheric  excitation  of  the 
plume,  is  shown  in  Figure  1-4.  The  viewer  angle,  plume  direction,  and  free 
stream  velocity  direction  are  all  completely  arbitrary.  (Section  3, 

Figures  3-1  and  3-2.  ) 

Only  one-half  the  plume  is  treated  by  a  single  FLAME  run.  The  remaining 
half  can  be  obtained  by  redefining  the  azimuthal  angles  as  shown  in  Section  3, 
Figure  3-2.  It  is  not  always  necessary  to  generate  both  plume  halves 
because  the  plume  radiance  is  axisymmetric  under  certain  conditions. 

The  radiative  models  used  tend  to  underpredict  the 

radiation.  For  atmospheric  excitation,  only  single  collisions  are  considered 
and  the  Boltzmann  distribution  of  excited  rotational  states  is  assumed.  Per- 
sistance  of  the  plume  radiance  due  to  pumping  by  thermal  (rather  than  pure 
kinetic)  collisions  with  the  atmosphere  is  not  treated.  For  particle  radiation, 
scattering  of  Earthshine  is  not  treated. 

Section  2  summarizes  the  theory  used.  Sections  3  and  4  detail  the  workings 
of  the  FLAME  code,  describing  the  equations  used. 
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Figure  1-3.  Isointensity  Lines  for  Look  at  45°  from  Plume  Axis  of  Solid  Propellant  Rocket 
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F  igure  1-4.  Isointensity  Lines  for  Broadside  Look  at  Particle  and  Excited  Water  Radiation 
from  Solid  Propellant  Rocket 
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Section  2 
THEORY 

The  major  elements  of  the  theory  used  as  the  basis  for  the  FLAME  code  is 
outlined  in  this  section.  The  approximate  model  of  the  far-field  multiphase 
plume  is  described.  For  two-phase  plumes,  closed-form  expressions  for 
particle  concentration  as  a  function  of  plume  coordinates  and  particle  size 
are  presented.  The  radiative  mechanism  for  the  particle  cloud  is  simply 
emission,  with  an  effective  minimum  temperature  defined  by  solar  illumin¬ 
ation  or  f  cattering  of  Earthshine;  it  need  not  be  elaborated  on  here.  Gas 
radiation  is  considered  to  be  caused  by  excitation  of  water  rotational  spec¬ 
trum  by  collision  with  the  atmosphere;  this  theory  is  outlined  in  subsection 
2.  3.  The  model  can  be  generalized  to  include  other  radiating  rotors. 

The  approaches  used  were  based  on  existing  theoretical  work.  No  attempt 
has  been  made  to  derive  new  theoretical  expressions  or  include  radiative 
mechanisms  not  yet  understood  in  the  current  FLAME.  The  program  is 
flexible  enough,  however,  to  allow  the  inclusion  of  other  mechanisms  at 
some  future  time. 

2.  1  PLUME  GAS  FLOW 

Closed-form  approximation  to  the  plume  gas  is  relied  upon  for  flow  applicable 
to  the  far-field  plume.  Atmospheric  mixing  and  the  resulting  plume  slow 
down  have  not  been  addressed.  Atmospheric  excitation  of  the  unaltered 
plume  is  considered. 

Closed-form  solutions  have  been  developed  for  the  density  profiles  as  a 
function  of  engine  operating  conditions.  It  appears  that  all  properties  of  the 
far-field  gaseous  plume  can  be  thus  treated.  The  computer  program  will 
respond  to  input  of  engine  thrust  (or  chamber  pressure),  chamber  tempera¬ 
ture,  expansion  ratio  (e),  Y  (or  exit  mach  number,  Mg),  and  nozzle  lip  angle 
by  generating  a  complete  history  of  the  rocket  plume. 
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To  predict  transition  to  free-molecular-flow,  a  relationship  was  derived  to 
exoress  the  Knudsen  number  (Kn)  variation  in  a  bipropellant  rocket  engine 
plume  as  a  function  of  the  polar  coordinates  (r,  0)  and  basic  engine  param¬ 
eters  F  (thrust),  Cp  (thrust  coefficient),  Tc  (combustion  temperature), 

M  (exit  Mach  number)  and  (ratio  of  specific  heats).  The  gas  viscosity, 

6  1/2 

[a,  was  assumed  to  vary  as  T  .  The  engine  chamber  pressure,  P  ,  could 
be  introduced  through  the  basic  thrust  equation: 

F  =  P  AC 
c  t  F 

where:'  A  =  Throat  Area 

These  derivations  used  Roberts’  approximation  (Reference  1)  for  the  plume 

2 

gas  density  which  indicated  that  the  density  was  inversely  proportional  to  r 

2 

and  proportional  to  cos  0  raised  to  the  power  Y(Y  -  1)  M^-2. 

The  derived  relationship  is  shown  in  Figure  2-1  as  a  function  of  Me>  Y  ,  and 
t  (engine  expansion  area  ratio).  This  information  may  be  employed  to  deter¬ 
mine  the  shape  of  the  surface  ( r,  0  )  required  to  achieve  a  given  Kn-  It  is 
noted  that  r  increases  with  increasing  Y  for  specified  values  of  T  ,  K  ,  F, 
and  6.  Indeed,  this  relationship  indicates  that  r  — 0  as  0  —  90  degrees,  i.  e.  , 
that  transition  occurs  at  the  nozzle  lip  where,  theoretically,  there  is  an 
infinite  pressure  gradient  as  the  flow  undergoes  a  Prandtl-Meyer  expansion. 

Figure  2-2  presents  the  variation  of  plume  density  on  a  spherical  cap  (r  = 
constant)  as  a  function  of  0  and  Mg  with  Y  =  1.4.  As  shown,  the  density  has 
been  normalized  by  its  value  on  the  engine  centerline,  i.  e.  ,  0=0.  When 

M  =2.  312,  the  density  distribution  is  linear  with  cos  0  and  becomes  increas- 

6  12 
ingly  dependent  on  cos  0  as  M  increases,  e.  g. ,  for  M  =  5.  0,  p  ~  (cos  0) 
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2.2  PARTICLES 

A  combined  nozzle  and  plume  flowfield  calculation  was  performed  for  an 
aluminized  solid  propellant  motor  to  obtain  correlations  of  particle  charac¬ 
teristics  in  the  near  field.  Closed-form  solutions,  analogous  to  those  obtained 
for  gas-phase-only  plumes  were  not  possible  due  to  the  nonequilibrium  nature 
of  the  expansion  involving  momentum  (velocity)  and  energy  (temperature) 
differences  between  the  phases. 


Correlations  were  obtained  which  related  particle  velocity,  temperature, 
density,  and  limiting  streamline  direction  to  the  seven  particle-size  groups 
considered.  An  expression  was  derived  which  related  the  particle  Knudsen 
number  to  plume  coordinates  and  motor  operating  conditions.  The  primary 
use  of  this  expression  was  to  determine  the  extent  of  the  particle  continuum 
region. 

Initial  estimates  of  particle  emissivity  were  made,  based  primarily  on  exist¬ 
ing  calculations  for  alumina.  The  details  of  these  analyses  are  presented  in 
later  subsections. 

2.2.1  Two-Phase  Flowfield  Particle  Correlations 

The  computer  programs  described  in  References  2  and  3  were  employed  to 
calculate  the  nozzle  and  plume  flowfields  for  an  aluminized  solid-propellant 
motor  with  the  following  characteristics: 

Thrust  =  15,000  lb^. 

Chamber  pressure  =  550  psia 
Expansion  area  ratio  =  30.  8 
Aluminum  weight  content  =16  percent 

A  total  of  seven  aluminum-oxide  particle  groups  were  utilized  with  particle 
radii,  pi,  ranging  from  0.  6  to  4.  7  pm.  The  corresponding  number  count 
versus  pj  was  obtained  from  Reference  4. 

The  computed  output  contained  particle  velocity,  temperature,  flow  directions, 
and  nondimensional  number  density  for  each  group  at  flowfield  mesh  points 
within  the  plume.  This  information  was  processed  to  extract  certain  particle 
asymptotic  characteristics,  namely;  velocity,  temperature,  and  flow  direc¬ 
tion  along  the  particle-limiting  streamlines.  (The  limiting  streamline  is 
defined  as  the  boundary  outside  of  which  a  given  particle  size  will  not  be 
present. )  It  was  found  that  these  streamlines  became  straight  lines  after 
the  particles  had  traversed  a  relatively  short  distance  due  to  the  rapid  decay 
of  interphase  momentum  and  energy  exchange. 
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The  values  of  particle  temperature,  T^,  flow  angle,  0  and  velocity,  V(P^,0), 
on  the  limiting  streamlines  are  shown  in  Figure  2-3  as  a  function  of  p^. 

It  is  noted  that  particles  larger  than  4jj.  are  still  undergoing  the  liquid-solid 
phase  transition  at  the  aluminum- oxide  melting  temperature  of  4,170  R  and 
that  these  particles  are  contained  within  a  23.  5-degree  halt-angle  cone.  Only 
forced  convection  heat  transfer  between  the  particles  and  the  gas  is  considered 
in  the  computer  programs  and,  therefore,  a  further  cooling  by  radiation  must 
be  considered  separately. 


The  plume  computer  output  was  employed  to  obtain  the  particle  velocities  for 
the  seven  discrete-size  groups  considered.  The  correlation  was  obtained  in 
the  form: 


v(p.,e  ) 


[a^  +  b^  exp  ( — c^  R/r^.)!  (cos  0)ni,  0  <  0  S  0^ 


(2-1) 


where:  R,0 

r 

'  t 
0 

Pi 


Polar  coordinates  in  plume 
Nozzle  throat  radius 

Particle  limiting  streamline  inclination  (Reference  1) 


Table  4-1  in  Section  4  shows  the  correlation  coefficients.  The  a.  coefficients 

i 

decrease  with  increasing  particle  diameter  which  indicates  that  the  larger 
particles  achieve  a  lower  limiting  velocity  as  R/r^.  »1.  In  addition,  the  m 
coefficients  were  all  negative  which  indicates  that  the  particle  velocities 
increase  off-axis  on  a  given  spherical  cap.  This  effect  is  due  to  the  higher 
particle  concentration  on-axis  which  has  caused  more  interphase  momentum 
transfer  from  the  gas  to  the  particles,  i.  e.  ,  on  a  per  particle  basis  the  gas 
is  less  about  to  "drag"  the  particle.  The  lower  off-axis  particle  density 
allows  the  gas  to  "drag"  the  particles  to  higher  velocities. 


The  particle  density  correlation  was  redone  specifically  for  the  off-axis 
dependency: 

dp .  r  n.  b. 

-- =  a  (^)  1  (cos  0)  \  0<6  <  0  (2-2) 

dPo  i  R  Pi 

The  correlation  coefficients  are  given  in  subsection  4.  5.  1.4. 
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Aside  for  the  largest  particle  size  considered,  the  correlation  showed  that 
b.  >0,  i.  e.  ,  the  particl«e  density  decreased  off-axis. 


The  two-phase  computer  results  were  analyzed  to  obtain  a  correlation  of  the 
gas  phase  mass  flux  on  the  engine  centerline  with  the  result  being: 


=  a 


-2. 7443 


(2-3) 


6=0 


where  d  =  chamber  density 
c 


The  above  result  indicates  that  the  mass  flux  is  decaying  more  rapidly  than 
inverse  square  with  distance  for  the  conditions  employed.  Additional  com¬ 
puter  runs  should  be  performed  for  lower  value  of  the  solids  loading  in  the 
exhaust  products  to  determine  if  the  gas  phase  mass  flux  will  approach  a 
^  dependency. 


Following  Reference  5,  the  off-axis  gas  phase  mass  flux  was  assumed  to  be: 


(vi)e 

(VDe  = 


(V  d 
where  c  = 


=  exp  [-X^(l-cosc0)^] 


(2-4) 


0 

90° 


6  +  (  V  -V  ) 

e  v  max  e' 


and  6g  =  Nozzle  lip  angle  at  exit 


v  -  Maximum  Prandtl-Meyer  angle 
max  /a 


vg  =  Prandtl-Meyer  angle  at  exit  Mach  number 


An  expression  was  derived  which  relates  the  particle  Knudsen  number,  Kn  , 

1 

to  engine  operating  conditions,  location,  and  particle  size  for  a  two- 
phase  plume.  This  expression  was  obtained  by  defining  Kn  as  the  ratio  of 

P{ 

the  relative  Mach  number  and  Reynolds  number,  i.  e.  ,  : 
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(2-5) 


Kn 


Pi 


Mr. 

1 


v-vDj|/yicYRT 

d|V-V  |2p.  /(JL 
Pi  Ki 


where: 

V  =  Gas  velocity 

V  =  Velocity  of  the  ith  particle 
Pi 


Y 

R 

T 

d 


Pi 


=  Conversion  factor  (=32.2  flbm/lbf] ).  (ft/sec2) 
=  Gas  ratio  of  specific  heats 
=  Gas  constant 
=  Gas  static  temperature 
=  Gas  density 

=  Radius  of  the  ith  particle 
=  Gas  viscosity 


Roberts'  solution  (Reference  1)  was  employed  to  relate  d  to  location  in  the 
plume  and  to  certain  engine  operating  parameters;  the  viscosity  was  assumed 
to  vary  as  \/f;  and  the  isentropic  relation  was  used  to  eliminate  T.  With 
these  substitutions  Equation  1  becomes: 

10 


Kn 


Pi 


12x46.  6x10 


T 

c  c 


P  Y  3/2  (Y- 1)  M  2 
c  e 


(1+  ¥  “a") 


1 

Y-l 


(r/r  V 
e 


P.  cos  CQ 
1 


wher  e:' 
r 

e 

R 


Me 
r,  e 
k 


=  Nozzle  exit  radius 


=  Universal  gas  constant  (=  1545  ft.  lbf/mole  °R) 
=  Combustion  temperature 
=  Chamber  pressure 
=  Nozzle  exit  mach  number 
=  Polar  coordinates  of  point  in  plume 
=  Hypersonic  similarity  parameter  (=Y(Y-l)m2) 


The  above  equation  indicates  that,  for  given  engine  operating  conditions  and 
plume  coordinates,  the  smaller  particles  will  achieve  the  higher  Knudsen 
number,  i.  e.  ,  they  will  "transition"  sooner.  Numerical  calculations  for  a 
lp  radius  particle  and  typical  engine  operating  conditions  indicates  that  Kn 
reaches  unity  very  close  to  the  nozzle  exit  plane.  By  comparison  a  gas- 
phase  plume  achieves  a  unity  Knudsen  number  only  at  large  distances  from 
the  exit  plane. 
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2.  2.  2  First  Estimate  of  Emissivity  of  Alumina  Particles 

The  Mie  theory  for  scattering  of  electromagnetic  waves  by  dielectric  spheres 

forms  the  basis  for  theoretical  determination  of  particle  emissivity.  The 

theory  actually  computes  the  scattering  and  absorption  efficiency  factors, 

Q  (a)  and  Q  (a),  defined  oy 
s  a 

Q  =  <r/o- 

g 

where  cr  is  the  geometrical  cross  section  and  cr  the  actual  cross  section. 

£> 

a  is  the  size  parameter  2Tra/\,  where  a  is  the  particle  radius.  The  values 
obtained  Jtor  Q  are  determined  by  the  values  of  the  real  and  imaginary  parts 
of  the  index  of  refraction  n,  related  by 

n  = 

Since  n  -  n(\,  T),  a  more  accurate  functional  dependence  of  Q  would  be 
Q=Q(\,  a,  T). 

Plass  (References  6  and  7)  has  computed  efficiency  factors  for  alumina 
spheres  of  0.  1  to  10p  for  a  range  of  wavelengths  0.  5  to  10p.  Values  of  n^ 
and  n?  used  in  this  wavelength  range  were  based  upon  measurements  made 
with  large  sapphire  crystals.  They  indicates  n^/n^  <<:** 

For  determination  of  the  emissivity,  e  ,  Plass  used  an  equation  originally 
derived  for  the  emissivity  of  a  semi-infinite  homogenous  slab,  with 

£=  Z-3\fj°~s 

Thus  the  following  values  of  «  are  obtained  from  the  following  table:. 


\ 

r 

€ 

2p 

1.  In 

0.  005 

5.  1 

0.  006 

9.  9 

0.  018 

5|x 

1.  1 

0.  054 

5.  1 

0.  089 

9.  9 

0.  115 

(Figures  for  \=  lOp.  are  not  immediately  available). 


It  is  to  be  noted  that  in  the  range  10<X<20p.,  values  of  n £  increase  by  more 
than  an  order  of  magnitude  (Brannon  and  Goldstein,  Reference  8)  so  e  would 
be  expected  to  be  larger  there,  aside  from  size  effects. 

The  results  of  Plass  have  been  questioned  by  Morizumi  and  Carpenter 
(Reference  9)  on  the  grounds  that  measurements  of  plume  emissivities  of 
rocket  exhaust  indicate  larger  values  for  particle  emissivities.  This  con¬ 
clusion  is  rationalized  on  the  grounds  that  absorbing  properties  of  rocket 
alumina  particles  are  higher  than  those  of  sapphire,  due  to  the  polycrystal¬ 
line  properties  of  the  former.  They  conclude  that  emissivity  of  alumina 
particles  in  rocket  exhausts  lies  between  0.  1  and  0.  3. 

It  is  concluded  that  conflicting  and  insufficient  data  available  at  this  time 
precludes  a  definitive  answer  for  the  emissivity  in  the  infrared.  For  pro¬ 
visional  values,  the  following  are  suggested,  using  linear  interpolation  and 
extrapolation  to  obtain  values  at  other  sizes  or  wavelengths. 


X 

r 

€ 

X 

r 

€ 

X 

r 

6 

5p 

lp. 

0.  05 

1  0|JL 

iM- 

0.  1 

2  Op 

iM- 

0.  2 

5 

0.  09 

5m- 

0.  15 

5 

0.  25 

10 

0.  12 

10  p 

0.  20 

10 

0.  30 

An  approximate  analytic  correlation  of  the  above  data  for  A^Oj  emissivity 
was  found.  The  correlation  is  shown  in  Figure  2-4,  and  is  adequate,  con¬ 
sidering  the  coarseness  of  the  available  data.  The  emissivity  correlation  is 
of  the  form 

«x(r)  =  (  X  +  2.  5)(.  00689  +  11.  32p  -  2220p2)  (2-6) 

for  X  in  microns  and  p  (particle  radius)  in  cm. 

2.  3  ATMOSPHERIC  EXCITATION 

Atmospheric  excitation  of  rocket-plume  species  by  collisional  high  relative 
velocities  has  been  identified  as  a  strong  source  of  plume  radiation.  A 
detailed  analysis  of  the  problem  was  presented  in  Reference  10  and  sum¬ 
marized  in  Reference  11.  The  unique  feature  of  this  analysis  is  the  inclusion 
of  the  time- dependent  radiative  decay  in  the  treatment  of  the  excited  species. 
Thus,  there  is  a  finite  time  required  before  the  plume  can  generate  radiation 
by  this  mechanism. 


Detailed  calculations  of  trajectories  of  the  colliding  species,  averaged  over 
initial  orientations,  were  done  using  Program  P2160.  Figure  2-5  shows  the 
resulting  rotational  temperatures  induced  as  a  function  of  impact  parameter 
and  Vco.  Also  shown  on  Figure  2-5  is  the  linear  variation  found  for  T^/Vco 
as  a  function  of  V«  (T  is  the  maximum  rotational  temperature  for  any 
impact  parameter). 


The  effective  collision  cross  section  (S)  to  use  with  is  found  by  requiring 

T  S  =  2^1”  T(b)  bdb 
m  ° 

Calculated  values  of  S  are  shown  ir  Figure  2-6.  For  the  current  version  of 
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FLAME,  S  will  be  assumed  to  be  constant  (=0.  7x1  0  cm  ). 

The  atmospheric  interaction  model  used  in  FLAME  follows  that  of 
Reference  10  namely: 

A.  For  the  wavelength  band  of  interest,  the  total  energy  (E^)  and 
radiance  (j)  per  molecule  are  calculated  as  a  function  of  T. 


(kfpj) 


C.  The  time  of  flight  to  a  position  in  the  plume  is  compared  to  the 
decay  time  so  that  an  average  radiance  per  excited  molecule  can 
be  found. 

D.  The  fraction  of  molecules  that  are  excited  are  calculated  assuming 
initially  cold  plume  species. 

E.  No  mechanism  currently  is  programmed  to  account  for  plume/ 
atmosphere  mixing,  but  plans  are  to  add  this  later.  Currently  if  a 
molecule  undergoes  multiple  collisions  in  the  atmosphere,  the 
resultant  radiance  is  found  by  simple  addition. 

F.  The  radiance  decay  of  an  excited  species  is  assumed  to  be  linear  in 
time. 

G.  Data  currently  in  program  is  applicable  to  a  single  species,  IT>0. 

Further  details  of  the  analysis  are  available  in  Reference  10. 


The  analysis  of  Reference  10  was  changed  in  one  way.  Since  the  radiative 
model  used  assumes  single  collisions,  and  the  plume  may  be  viewed  from 
long  distance,  an  attenuation  factor  (A^)  must  be  introduced  to  account  for 
depletion  of  the  supply  of  uncollided  molecules  as  the  gas  expands  far  from 
the  nozzle.  Assuming  each  collision  removes  that  molecule  from  further 
consideration,  the  reduction  of  the  number  of  available  molecules  from  NQ 
(initial  value)  to  N  (at  point  r  in  plume)  is  given  by 


exp  (-At) 


where 


(2-7) 


r 

At  -  s  /  nwds, 

o 

where  S  is  the  collision  cross  section  and  n^  is  the  atmospheric  density. 


A  correlation  was  obtained  for  the  spectra)  linear  absorption  coefficient  of 
the  H^O  rotational  spectrum.  The  results  are  that 


=  log-  ^  (y) 


-1 

in  cm 


(2-8) 


wher  e 


2 


-3. 208 


279.  96 

yo  \fT  4-  50.  031 
A  -  [1.030  -  -*°5,582  ]  xlO’4 

JT  -2.896 

w  =  10,000/X.,  X.  in  pm 


The  data  correlated  are  shown  in  Figure  2-7,  where  the  results  of  the 
correlation  are  shown  for  temperatures  of  600°K  and  1,800°K. 
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Section  3 
FLAME  CODE 


The  following  subsections  deal  with  the  logic  underlying  the  FLAME  code  as 
it  is  now  configured.  Subsection  3.  1  describes  the  logic  that  went  into  the 
choice  of  coordinate  systems,  and  displays  the  geometric  model  used.  Sub¬ 
section  3.2  presents  the  flow  charts  used  as  the  basis  for  FLAME. 

3.  1  GEOMETRY  CONSIDERATIONS 

The  rocket  plume  is  an  axisymmetric  body  about  the  z-axis,  as  illustrated  in 
Figure  3-1.  Integration  through  the  plume  is  carried  out  through  the  conical 
volume  in  the  first  and  second  quadrants  defined  by  the  polar  coordinates  0m 
and  r  .  The  viewing  angle  is  specified  completely  by  the  angle  a  with  respect 
to  the  x-axis  in  the  xz  plane. 

We  want  to  integrate  emittance  along  families  of  lines  of  sight,  such  as  s, 

/ E(r,  0  )  ds 

We  can  express  ds  in  terms  of  either  cartesian  or  spherical  coordinates, 

ds  =  dx 
or 

ds  =  dr 

Since  emittance  is  a  function  of  r,  0  along,  equation  (3-2)  is  more  inviting  to 
use,  if  it  is  practical.  The  following  analysis  shows,  however,  that  cartesian 
coordinates  will  probably  result  in  less  computation  time,  and  should  be  used. 


The  conical  segment  in  Figure  3-1  has  equation 

2  ,  2  2.  2C  ,, 

x  +  y  =  z  tan  0j  (3-3) 

The  family  of  lines  along  direction  of  observation,  a .  have  equations 

y  =  yx  (3-4) 

x  =  Xj  +  s  cos  0  (3-5) 

z  =  z^  +  s  sin  a  (3-6) 
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>♦(£)*♦(£) 2 

!+r2(|.f)2+r2sin24.(|£)2  (3.2) 


Adding  Equations  2-10  and  3-11 


2  2  2  2  2 
Yj  +  cota  (rcosO-Zj)  +  =  (rsinG)  =  r  -  (rcosG) 


cot  a 

r  cosG  = 

F 

-  cot  O' 

zi  +  xi  = 

A 

2 

2 

-Yl  + 

r  =  B 

Equation  3-12  becomes 


(3-12) 


?  2  F 

F  +  2AF  +  A  -  B  +  -£-~ ^  =  0 

cot  o 


F2(l+tan2a)  +  2AF  +  A2  -  B  =  0;  F2/cos2o  +  2AF  +  A2-B  =  0 


Ther  efore 


(  1  \ 

1 4A2-4A2  | 

2  ) 
(cos  o/ 

1  /  cos 


2  /  2  2 
•  2A  cos  a  ±  2  coso  \  B-A  sin  a 


Therefore 


/  2  2 

sin  a  r  cosG  =  -2A  cos  a  ±  2  \  B-A  sin  o 


cos  6  = 


2Acot  a  _j_  2, 
r  r 


/-|--A2  su(r) 
/sin  a 


„  -l/2Acota  ,2  /  B  <2 

e  =  cos - —*7  \  TT'a 

\  V  sin  a 
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From  consideration  of  the  case  a  -  90  degrees  and  noting  that  cos6>  0,  only 
+  root  is  valid. 


d0  - 1  du 


dr 


\R 


dr 


or 


d0  _  -1 

dr 


vi-u2 


/  2\ 

-i/2 

_  u.  -f. 

I 

'  B  -A  > 

2r 

r 

r 

i.2 

•  2„ 

\sm  a  / 

sin  o . 

also 


sin4>  = 


1 


r  sin0 


v/(l-u2) 


=  V 


(r) 


d<|> 
dr  " 


1  _  dv 

f 


where 


dv  _  ~yl  ,  yl  U  du 

dr  "  r^  /T72  r(l-u2)3/2  dr 


(3-13) 


From  Equation  3-13 

dG  _  uw 

d?  ~  2.1/2 

r  (1-u  ) 

where 


w(r)  = 


1  - 


2r 


.  2 

usin  a 


^  /  s  in2  n  -  A2 


then 


dv 

dr 


Y,  u  w 


2m  2, 
r  (1-u  ) 


372 
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; 


du 

<Jr 


uw 

r 


or  finally 


ds  =  dr  1  +  r 


-2  _J_  yI  /  +  ArV  +  ^  yl  ul  2 

1-v2  r4(l-u2)  \  1-u2/  r2(l-u2)  r2 


(3-14) 


From  Equations  3-3,  3-4,  3-5,  and  3-6  the  limits  of  integration,  r,  and  r^ 


can  be  found. 


We  have 


2  2  ,  2  2 
x  +  y  +  z  =  r 


(3-15) 


2,2  2t  2  22 
x  +  y  =  z  tan  o  ,  =  r  -z 
7  °  m 


2  2  ,  2  , 
r  "  z  (1  +  tan  0  ) 

m 


Ther  efor  e 


z  =  rcos0 


Also 


Y  =  Yi 


X-X|  z-z^ 
coso  sinn 


X  =  coto  (z-z^)  +  Xj 


(3-16) 
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Substituting  x,  y,  and  z  into  equation  of  the  cone.  Equation  3-3  we  have  a 
quadratic  equation  for  r,  namely 


2  2  2  2  2 
(cot  a  cos  Gj  +  cos  0  ^  -l)r  +  (2x^  cotacosG^  -2  cot  o  cos  0  j)r 

+  (cot2o  z2^  -2x^coto'z^  +  x^2  +  yj2)  =  0  (3-17) 

The  solution  gives  r  ^  and  r^,  as  function  of  z^,  x^,  and  y^ 


3.1,2  Cartesian  Coordinates 
From  Equation  3-16 

z  =  tan  a(x-x^)  +  z^  (3-18) 


or 


dz 

dx 


t?n  u 


Therefore 


ds 


dx 


1  +  tan  a 


dx 

cos  a 


(3-19) 


or 


ds 


dz 

sina 


(3-20) 


We  obtain  r  and  0  from  Equation  3-18 

r2  =  z2  +  y^2  +  jcot  a  (z-z^)  +  x^j  (3-21) 

cos  6  =  z/r  (3-22) 

It  is  apparent  that  it  is  easier  to  use  cartesian  coordinates,  and 
Equations  3-19  or  3-20  together  with  Equations  3-18,  3-21,  and  3-22  to 
integrate  through  the  plume,  than  to  use  the  single,  complex  Equation  3-14. 


To  avoid  singularities  at  a  =  0  degree  or  90  degrees,  x  will  be  used  as  the 
independent  variable  for  0  <  a  <  45  degrees  and  z  for  the  independent  variable 
for  45  degrees  <  a  <  90  degrees.  Since  the  plume  is  transparent,  negative 
values  of  a  are  not  needed,  the  integrals  being  equal  to  values  for  90  degrees 
+  a . 

For  cartesian  coordinates,  the  intercept  equations  are  more  complex  than 
for  spherical  coordinates;  substituting 


y  = 


X  -  COtO  (z-Zj)  +  Xj 


into  equation  of  spherical  cap 


2,2,2  2 

x  rv  +  =  r 

m 


and  conical  sides 


(3-23) 


(3-24) 


(3-25) 


jj|fl 


2  2  2  2 
x  +  y  =  z  tan  0 


(3-26) 


yields  the  following  equations  for  the  resulting  four  intercept  points  of  a 
straight  line  of  sight: 


For  0  <  a<  45  degrees  (choose  x^  =  0) 

(1  +  tan ^o)x^  +  (2  z^  tana)  x  +  -rm^  =  0 

(cot^0 m-tan^»  )x^-(2zj  tana)  x  +  y^cot^O^  ~z\  ~  ® 


and  for  45  degrees <  a  <  90  degrees  (choose  Zj  =  0) 

(l-cot^a)z^  +  (2x,cota)z  +  y.^  +  x/-r  ^  =  0 
1  1  1  m 

(tan^0m~cot^a  )z^-(2xjCota)z  -  ^y.^  +  =  0 


(3-27) 


(3-28) 


(3-29) 


(3-30) 


For  each  range  ofa,  solution  of  the  two  quadratic  equations  will  yield  four 
values  for  x  or  z.  The  correct  limits  of  integration  are  found  by  selecting 
those  pairs  of  x  or  z  for  which  the  corresponding  values  of  r  and  0  (by 
Equations  3-18,  3-21,  and  3-22)  are  on  the  real  plume  boundary,  namely 

r  and  8  <  0 
m  m 


0  and  r  <  r 
m  '  m 


3.1.3  Integration  Schemes  and  Step  Sizes 
Most  of  the  integrals  of  the  form 


I  =  f  f(x)< 


(3-31) 


are  evaluated  using  Gaussian  eight-point  quadrature.  Defining 


y  =  (2x  '  x2  -  *  xj ) 


(3-32) 


we  have 


I  =  X2  X1  /  f(x)dy 


(3-33) 


Then 


o 

x2’xl  V 

i  2L  w-  f<x-) 

2  f— '  .  i  l 


(3-34) 


where 


=  0. 18343  46424  95650 
=  0. 52553  24099  16329 


=  0.79666  64774  13627 


=  0.  96028  98564  97536 


and 


w.  =  w0  =  0.  36268  37833  78362 
1  o 

w2  =  w?  =  0.31370  66458  77887 

w3  =  w6  =  0.22238  10344  53374 

w  =  w  =  0.  10122  85362  90376 
4  5 

The  values  of  are  related  to  through  Equation  3-32. 

For  integrating  through  the  plume,-  Simpson's  rule  is  used 
x+2Ax 

f(y)dy  =  (f(x)  +  4f(x+Ax)  +  f(x+2Ax))  (3-35) 

x 

The  step  size,  Ax,  is  varied  after  each  successive  integration  from  x  to 
x  +  2Ax.  A  geometi  ic  variation  of  Ax  assures  that  the  smallest  step  size 
occurs  as  the  integration  along  the  line  of  sight  passes  through  the  plume 
center  (x  =  0  or  z  =  0  planes). 

Let  the  total  maximum  distance  along  the  line  of  sight  be  x.p.  Starting  at 
the  middle  of  the  plume  and  integrating  out,  a  distance  x^,/2  is  covered  in 
about  n/2  steps  in  is  a  program  input).  If  each  Ax  step  is  r  times  bigger  than 
the  previous,  the  rule  for  geometric  progression  gives 


Ax 


min 


-1 

1 


(?-36) 


where  Ax  .  is  the  initial, 
mm 

increase  in  Ax  each  step,  r 


smallest  step  size. 
-  1.  I 


Choosing  a  10-percent 


or 


Ax  . 


min 


20(1.  ln/2-l) 


(3-37) 
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Starting  at  position  x  =  -  |  xq|,  a  large  step  size  i.s  used  which  decreases  until 
x  =  0  and  then  increases  again.  The  initial  value  of  Ax  is  found  from  the 
formula  for  the  sum  in  terms  of  the  last  term  (Ax^n) 


x 

o 


(3-38) 


or 


Ax.  =  (0.  1  x  +  Ax  .  )/l.  1 
m  c  nun 


(3-39) 


A  similar  stepping  scheme  is  used  to  obtain  a  higher  density  of  lines  of  sight 
near  the  rocket  exit  plane  or  plume  centerline  than  used  in  the  rarified  plume. 


3.1.4  Relative  Velocity 

The  rocket-body  velocity  must  be  specified  in  the  coordinate  system  of 
Figure  3-1  so  that  relative  velocity  of  plume  gases  and  atmosphere  can  be 
calculated.  Figure  3-2  shows  the  position  of  the  vehicle  velocity  vector 
V^,  and  the  definition  of  straight  up,  u. 

Integration  is  performed  only  over  a  plume  where  y  >  0.  To  obtain  other 
half  of  plume,  and  should  be  increased  180°. 

Relative  velocity  is  just  the  vector  sum  of  Vpr  and  VR  where  r  is  the  unit 
vector  along  the  polar  coordinate  r  to  a  point  in  the  plume  being  considered. 
Altitude  at  this  same  point  is  missile  altitude  hg  plus  projection  of  7"  onto  u. 

h  =  hQ  +  7-  u  (3-40) 

3.  2  FLOW  CHARTS 

This  section  presents  the  logic  flow  of  the  FLAME  code.  Input  data  are 
listed  in  Tables  3-1  through  3-4.  The  Figure  3-3  is  a  simplified  diagram  of 
the  overall  code.  The  charts  (Figures  3-4  through  3-9)  expand  upon  certain 
regions  of  the  computer. 
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Table  3-1 

LISTING  OF  INPUT  DATA 


Read  in 


Real 

or 

FORTRAN  Integer 


t. 

time  after  ignition  (sec) 

TIG 

R 

r  , 

0 

nozzle  exit  radius  (cm) 

RE 

R 

Y, 

ratio  of  specific  heats  (>1) 

GAMMA 

R 

M  , 

0 

exit  mach  number 

ME 

R 

V 

nozzle  lip  angle  (deg) 

LIP 

R 

de(t). 

exit  plane  density  as  a  table  versus 

PDE 

R 

time  after  ignition  (g/cm^  vs  sec) 

TID 

R 

!nD. 

night/day  flag  (0  or  1) 

IND 

I 

tSL> 

solid /liquid  flag 

0  or  -2  (gas);  1  or  -1  (solid) 

ISL 

I 

Xmin>1 

WLMIN 

R 

\naxd 

band  limits  (pm) 

WLMAX 

R 

n, 

approximate  integration  steps  through 
layer 

INT 

I 

m. 

approximate  number  of  lines  of  sight 
used  in  one  direction 

LINES 

I 

Tc 

chamber  temperature  (°K) 

TC 

R 

xrrv 

radiating  gaseous  species  mole  fraction* 

XM 

R 

TEj 

Earth  effective  temperature  (°K) 

TE 

R 

look  angle  (0  <  o  <  90  degrees) 

AL 

R 

h0’ 

altitude  of  engine  (ft)  (from  0.  3  Mft  to 

3  Mft) 

HO 

R 

VR, 

rocket  velocity  (fps) 

VR 

R 

^u 

polar  angle  of  vertical  direction  (deg) 

PHIU 

R 

^u 

azimuthal  angle  of  vertical  direction  (deg) 

XIU 

R 

<^>R, 

polar  angle  of  rocket-body  velocity  (deg) 

PHIR 

R 

^R> 

azimuthal  angle  of  rocket-body  velocity 
(deg) 

XIR 

R 

END, 

flag  signifying  end  of  run  (0  or  1) 

END 

I 

*FLAME  assumes  an  average  molecular  weight  (WMOL)  of  17  for  the 
plume  species  (approximately  correct  for  pure  gas  plumes,  or  the 
gaseous  constituents  of  solid- propellant  plumes).  When  gaseous  radia¬ 
tion  from  a  solid  propellant  plume  is  calculated  (  |IgjJ  =  1  and  xm  ^  0)» 
the  molecular  weight  should  be  corrected  to  include  solids.  This  is 
done  by  modifying  xm,  since  only  the  ratio  XM/WMOL  is  used  by 
FLAME.  For  a  typical  solid-propellant  plume,  WMOL  =  26,  so  XM 
should  be  reduced  by  the  ratio  (17/26). 
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Table  3-2 


PROGRAM  REGION  1.  INITIAL  CALCULATIONS- PARTICLES 


The  following  data  will  be  calculated  for  each  particle 
radius,  p^,  corresponding  to  the  ith  particle  size: 


EPT, 

TCOOL, 

BPASS, 

ETL 

TMIN, 

DELT, 


table  of  total  emissivity  versus  temperature 

table  of  particle  temperature  versus  time  after 
solidification 

black  body  intensity  in  band  of  interest 

table  of  average  emissivity  in  region  from^min 
to  ^rnax  versus  temperature 

equilibrium  temperature  of  particle 

time  required  to  solidify  droplet 


Table  3-3 

PROGRAM  REGION  2.  INITIAL  CALCULATIONS-GAS 

Calculate  limiting  plume- gas  velocity,  VP. 

For  the  band  width  specified,  calculate 

ET,  table  of  total  available  rotational  energy  as 

functior  of  temperature 

JM,  table  of  radiative  rate  (radiance  per  molecule) 

as  function  of  temperature 

JT,  table  of  radiative  rate  as  a  functior?  of  time 


Table  3-4 

PROGRAM  REGION  3.  SELECT  X!  OR  zx  for  o  >  90  -  6m 

Test  IA  to  determine  magnitude  of  a.  For  LA  =  0,  define  zj. 

For  IA  =  1,  define  x-^.  Both  xj  and  zj  are  labeled  XI  in  the  program. 
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CALL  STAB  AND  CHECK  (YIELDS  INITIAL  AND  FINAL  VALUES  FOR 
INTEGRATION  THROUGH  PLUME) 


DEFINE  ANGULAR  QUANTITIES.  ONCE 
ONLY 


INOUT  *  0 


NEXT  VALUE 


CALCULATE  X  OR  Z  INTERCEPT  VALUES 


CALL  CHECK 

CHECK  TO  SEE  IF  IN  PLUME 


NOT  IN 
PLUME 


IN 

PLUME 


IS  THIS  FOURTH  VALUE? 


IS  THIS  SECOND 
INTERCEPT? 


INOUT  «  ^ 


DEFINE  PMIN  AND 
PMAX 


RETURN 

RETURN 

Figure  3-4.  Program  Region  4  Flow 

CALL  FUN  (EVALUATES  FUNCTION  BEING  INTEGRATED  THROUGH  PLUME) 


^  *0 

FIND  PARTICLE  TEMPERATURE 


RETURN 


Figure  3-5.  Program  Region  5  Flow 


CALL  EMIT  (R,  WL) 


IS  THIS  FIRST  TIME  EMIT  WAS  CALLED? 


HAS  p,  CHANGED  FROM  PREVIOUS  CALCULATION? 


CALCULATE  p-DEPENDENT  PART  OF  CORRELATION 


RETURN 


Figure  3-7.  Particle  Emissivity  Routine  Flow  Chart 


CALL  ABCO  (WL,  TE) _ 

HAS  TEMPERATURE  CHANGED  FROM  PREVIOUS 
CALCULATION? 


CALCULA  FE  TEMPERATURE 
DEPENDENT  PART  OF  k* 


CONVERT  k  TO  w 


CALCULATE  k* 


RETURN 


Figure  3-8.  Absorption  Coefficient  Routine  Flow  Chert 
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CR156 


CALL  CONTOUR  (KP> 


CHOOSE  ODD  NUMBER  OF  POINTS 
TO  BE  PLOTTED 


PLOT  GRID,  IF  NOT  YET  DONE 


FIT  2  DIMENSIONAL  SURFACE  TO  J<*.  V> 


FIND  LIMITS  OF  INTENSITY  IN  3  X  3  MATRIX 


DETERMINE  IF  A  GIVEN  VALUE  OF  ISOINTENSITY  TO 
BE  PLOTTED  IS  IN  THE  LIMITS  OF  THIS  MATRIX 


DETERMINE  COORDINATES  OF  20  POINTS 
OF  THIS  ISOINTENSITY  LINE  WITHIN  THE 
MATRIX  LIMITS  USING  CURVE  FIT 


NEXT  VALUE _ 

OF  ISOINTENSITY  LINE 


PLOT  LINE 


LAST  VALUE;  NEXT  MATRIX 


LAST  VALUE,  LAST  GROUP  OF  3  X  3  MATRICES 


PRINT  COORDINATES  FOR  LARGEST  ISOINTENSITY 
LINE  IN  ThiS  COLUMN,  IF  NOT  ALREADY  DONE 


INITIALIZE  FOR  NEXT  GROUP 
OF  MATRICES 


RETURN 


Figure  3-9.  Contour  Plottar  Flow  Chart 
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Section  4 

EQUATIONS  TO  BE  SOLVED 


The  following  subsections  relate  to  the  program  regions  shown  on  the  master 
flow  chart.  Figure  3-3. 

4.  1  PROGRAM  SECTION  1.  INITIAL  CALCULATIONS— PARTICLES 


4.  1.  1  EPS,  Particle  Emissivity  Versus  Temperature 

A 


tix  (T)  dx/ff  T- 


where 


R°  (T)  .  && 


he 

-1 

r* 

'  C2 

eikt  -  1 

cl 

e"T  -  1 

'  X5 

_ 

•  * 

1  -1 


c  =  velocity  of  light 
a  =  Stefan- Boltzmann  constant 

=  5.  672  x  10~^  j/cm^  °K^  sec 
eiA  =  correlation  eq  (2-6) 

Cj  =  3.  742  x  10^  w  (|im)Vm^ 
c2  =  14390  pm  °K 

See  Subsection  3.1.3  for  description  of  Gaussian  Integration. 


(4-1) 


(4-2) 


4.  1.2  TCOOL,  Particle  Temperature  Versus  Time  After  Solidification 


c  n. 

At;  (T)  = 


•  i  r  1  dT 

*  /  t.  T4 
J  T/  1 


(4-3) 


where 

Cp  =  particle  specific  heat 
Tf  =  fusion  temperature 


4.  1.  3  BPASS  Black  Body  Intensity  in  Band  of  Interest 
,  r  ^max 

Ba°(T)  =  i  j  R°(T)dK 

'  ^min 

4.  1.  4  EPSI,  Average  Particle  Emissivity 


‘i.AX  (T>  =  7 


‘max 


c  Kl  (T)  dX/B°x  (T) 


mm 


4.  1.  5  TMIN,  Particle  Equilibrium  Temperature 
at  day. 


4  i1  i  4 

(.  T  .  .  =  -r —  +  -r-  T_,  =  K 
i  mm,i  2  a  2  E 

where 

cu  =  solar  absorptivity  [~€.  (6000°)] 

a'.  absorptivity  to  Earthshine  [~e^  (T^,)] 

2 

T  =  solar  constant  =  0.  1353  w/cm 
Tj,  =  effective  Earth  temperature 


at  night.  Equation  4-6  is  used  with 

r  =  0 


Newton-Rapheson  solution  is  of  form 


T. 

J 


1_ 

3 


2T.  ,  + 
3-1 


K 


T3  ,  e.  (T.  .) 
3-1  i  3-1  . 


4.  1.  6  DELT,  Time  to  Solidify 


At. 

i,s 


p.nU 
3  c.  <r  T, 


where 

m  =  particle  density 
\  =  particle  heat  of  fusion 


(4-4) 


(4-5) 


(4-6) 


(4-7) 


(4-8) 
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4.2  PROGRAM  SECTION  2.  INITIAL  CALCULATIONS— GAS 


4.  2.  1  Plume  Gas  Limiting  Velocity 


/  2  \1/2 

V  = 

a  (  ,  ) 

(4-9) 

p 

°v  -  v 

where 

a  = 

\/  YRT 

(4-10) 

o 

c 

4.2.2  ET 

,  Rotational  Energy  Versus  Temperature 

et  b 

hcB  F  1/2  -xi  1/2  -X2 

+  Vj  (erf  ^ 

-  erf  '*/xj) j 

-3/2  X1  C  "  X2  C 

CO 

(4-11) 

where 

B 

=  molecular  rotational  constant 

a 

=  hcB/kT 

(4-12) 

Xl,2 

.  2 
=  J1>2  a 

(4-12) 

T 

1 

(4-14) 

1,2 

"  2BA  . 

mm,  max 

he 

2 

=  0. 9928  w  sec  cm 

he 

k 

=  1.4387  cm  °K 

4.2.3  JM 

,  Radiative  Rate  in  Band 

i  max  n 

(4-15) 

j(T)  = 

-  i  /  ^  <T>  <  iT) 

where  k^(T)  =  the  spectral  linear  absorption  coefficient  per  molecule  from 
Equation  2-8,  multiplied  by  kT/p  to  change  from  per  atm. 


4.2.4  JT,  Radiative  Rate  Versus  Time 
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where  T  is  found  for  value  of  E^,  from  results  of  S\ib  section  4.  2.  2  above,  and 
then  j  is  formed  from  results  of  Subsection  4.  2.  3. 


4.  2.  5  Closed  Form  tor  Radiance 

A  least- square  fit  is  made  to  the  resulting  data  of  Subsection  4,  2.  4  of  the 
form 

j(t)  =  A  ( tr  -  t)  (4-17) 


Substituting  Equation  4-17  into  4-3  6  yields 

T  2  ET 

r  ZttA 


(4-18) 


which  relates  to  T  ,  the  maximum  temperature  excited  by  atmo spheric 
colli  sion. 


4.  3  PROGRAM  SECTION  3.  SELECTION  OF  x  OR  z^ 

These  are  the  initial  values  of  or  z^  to  be  used  to  define  lines  of  sight  just 
grazing  plume. 

For  a  <  45  degrees,  we  have  the  following: 


6  <90  and  90-6  >  a 

m  m 


2  = 

l 

0 

(4- 

•19a) 

6  >90  and 

e  - 

90  <  90  -  a 

m 

m 

Z1  = 

-r  /cos  a 
m 

<4- 

•  19b) 

Otherwise 

II 

r 

m 

cos 

(a  +  0  )/cos  a 
m 

14- 

•19c) 

For  a  >  45  degrees,  we  have  the  following: 

0  <90  degrees  and  a  >  6 

m  &  m 

x.  =  -r  cos  (a  -  0  )/sin  a  '  (4-20a) 

i  m  m 


Otherwise 


In  addition,  if  there  is  no  gaseous  radiation  (XM  =  0),  the  plume  angular 
extent  is  taken  equal  to  the  flow  angle  of  the  limiting  streamline  for  the 
smallest  particle  (33.  52  degrees). 


4.  4  PROGRAM  SECTION  4.  CALL  STAB --INTERCEPT  VALUES 
For  0  ^  a  <  45  degrees 

2  ?  2  2  ? 


m 


and  for  45  degrees  <  a  <  90  degrees 

2  2  2  2 
(1  -  cot  a)  z  +  (2  x^  cot  a)  z  +  y^  +  x^ 


r  “  =  0 
m 

(3-27, 
Section  3) 

cot2  6  -  z  2  =  0 

(3-28, 

m  1 

Section  3) 

-r  2  =  0 

(3-29, 

in 

Section  3) 

:  +  Xj2)  =  0 

(3-30, 
Section  3) 

4.  5  PROGRAM  SECTION  5.  CALL  FUN— RADIANCE  CALCULATIONS 


z  =  tan  a  (x  -  x^)  +  z^ 


=  yjz2  +  y2  +  [cot  a  (z  -  Zj)  + 


0  =  cos  *  (z/r) 
4.  5.  1  Solid  Particles 


(3-18, 
Section  3) 

(4-21) 

(4-22) 


4.  5.  1.  1  Calculate  Throat  Radius 


rt  =  re/'jf? 


where 


Y+l 


Y+l 


‘  -  Jr ^  M°)  2<V‘1) 


(4-23) 
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4.  5.  1.2  Particle  Limiting  Velocity 

V(Pi  6)  =  fai  +  exP  ^-ci  r/r>l  c°sni  0  (4-24) 

for  0  <  9  <  0^  (limiting  particle  streamline  inclination)  where  the  coeffi¬ 
cients  are  given  in  Table  4-1. 


4.  5  .  1  ,  3  Particle  Time  of  Flight 


t  =  r/V 


(pvo) 


(4-25) 


4.  5.  1.4  Concentration 


d.  /r  \  ei 

Density  ratio  =  --p  =  f  (p^,  r,  0)  =  g.,  (— I  (cos  6) 


li 


(4-26) 


for  0  <  0  <  0^  and  the  coefficients  are  given  in  Table  4-2.  Concentration 
is  related  to  density  by 


N 


o  =  V(!K3  d'  <XML)-  d'  =  alumina  density 

=  4.  0045  g/cm^ 

XML  =  alumina  mass  fraction 


4.  5.  1.  5  Particle  Temperature 

If  t^  <  AL  g  (Equation  4-8), 

T  =  Tf 

If  tf  >  At.  , 
f  i,s’ 

T  =  maximum  of 

a.  TCOOL  (t,  -  At.  ),  Equation  4-3 

X  1  j  s 

or  b.  TMIN,  Equations  4-6  or  4-7 
4.  5.  1.  6  Evaluate  Integrand 

Ji  =  No  (t'  '  V  f(pi'  r’  0)  *  Pi  ei,A\  (T)  BA\  (T)  (4-27) 

(See  Equations  4-26,  4-4,  and  4-5) 
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Table  4-1 


PARTICLE  LIMITING  VELOCITY  COEFFICIENTS 


Pi(pm) 

0 . 

1 

a.(cm/ sec) 

b.  (cm/sec) 

l 

c. 

l 

n. 

l 

0.  690 

33. 52° 

30. 776  x  104 

-4.  6692  x  104 

0.01874 

-0.9048 

1.  330 

28.  22 

29.701 

-3. 8094 

0.01930 

-1.  3494 

2.  000 

23.  53 

28.932 

-3. 2891 

0.01976 

-1.  5419 

2.  660 

19.  76 

28.248 

-3.0367 

0. 02012 

-1.4886 

3.  332 

16.  81 

27. 520 

-2.9355 

0.  02036 

-1. 3175 

4.  000 

14.  75 

26.  764 

-2.9337 

0. 02064 

-1.0308 

4.  666 

13.  56 

26.007 

-2.9745 

0. 02092 

-0. 6425 

Table  4-2 

CONCENTRATION  COEFFICIENTS 

Pi(|J-m) 

gi 

e. 

l 

f. 

l 

0.  690 

0. 2201 

2.  33 

2.  84 

1.  330 

0.  4503 

2.24 

5.  64 

2.000 

0.  5395 

2.  17 

6.  75 

2.  660 

1.0308 

2.  12 

7.  30 

3.  332 

1.3339 

2.08 

6.95 

4.000 

0. 8870 

2.  (.  j 

6.30 

4.  666 

0.  5396 

2.04 

-2.  05 

4.  5.  2  Excited  Gases 

4.  5.  2.  1  Time  of  Flight 

t^  =  r/V  ,  Equation  4-9 

4.  5.  2.  2  Gas  Concentration 

k/re\2  .  ,  flUk-2 

n  =  n  — )  [cos  (a0)l 

p  p,e  2  \  r  / 


(4-28) 


(4-29) 
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where 

k 


Y  (Y  -  1) 


(4-30) 


a 


V 

max 


90  degrees 

6  +  ( v  -  v  ) 

e  max  e 


Y  +  1 

Y  -  1 


-  1 


90  degrees 


(4-31) 


14-32) 


1  Vttt  <M«-1)  -lan'1  Ar^-1  ,4-33) 

and  n  =  d  N  /M 
p,e  e  o 

Nq  =  Avogadro' s  number 
M  =  mean  molecular  weight 
and  remaining  parameters  are  input  data 


4.  5.  2.  3  Relative  Velocity 


x  =  cot  a  (z  -  z^)  + 


V« 


where 

VI 

V2 

V3 


1/2 

(4-34) 

v r  sin  «j>R  cos  vjjR 
v r  sin  <t>R  sin  i|jr 

VR  COS  +R 


4.  5.  2.  4  Excitation  Temperature 

T  =  maximum  of 
m 

A.  T  of  ambient  atmosphere  (See  Equation  4-40) 

B.  Tm  =  9.61194  x  10"9  V.2  +  1. 07636  x  10"3  Vm  (4-35) 

(Vco  in  cm/sec) 
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4.  6  PROGRAM  SECTION  6.  CALL  SIMP— INTEGRATION  THROUGH  PLUME 
Simpson's  rule 

f(x)  dx  =  ~  [f(s,  +  4  f^s  +  y)  +  f ( s  +  ds)]  (4-42) 

Step  size  is  controlled  by  a  geometric  progression,  so  that  ea  h  step  is 
10  -  larger  or  smaller  than  the  previous,  and  minimum  step  size  occurs  at 
z  =  0  or  x  =  0  planes. 


4.  7  CURVE  FITTING  ROUTINE,  CF 


An  auxiliary  routine  was  written  to  obtain  a  curve  fit  to  the  intensity  data  to 
facilitate,  plotting.  The  curve  fit  for  intensity  is  of  the  form 


•z 

i=l 


a^  (x,  y) 


where 


♦j 


2  2 

x  ,  x,  y  ,  y,  xy. 


1 


The  least  square  fit  condition  is  (for  nine  points) 


da 


Z  •  Z  aj  *j  <3V 

i=l  3=1 

The  resulting  matrix  of  equations  is 


k  / 


=  0 


(C)  (a)  =  (I) 
where 
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and  each  term  in  the  C  and  I  matrices  is  summed  over  all  points,  e.g. 

9 


i=  1 


This  technique  of  extrapolation  is  similar  to  that  used  by  Dailey  (Reference  12) 
Nine  points  (3x3  square)  are  selected  and  the  curve  fit  made.  Coordinates 
of  the  isointensity  lines  are  then  deduced  from  the  curve  fit  by  solving  the 
quadratic  equation  for  y,  given  I  and  x.  These  lines  are  plotted,  and  then 
the  next  group  of  nine  points  are  selected.  The  procedure  is  repeated  until 
the  plume  is  covered. 
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Section  5 

SAMPLE  PROBLEM  AND  OPERATING  INSTRUCTIONS 

A  list  of  all  input  quantities  is  given  in  Table  3-l*in  Subsection  3.  2,  with 
their  appropriate  units.  Additional  description  of  the  input  parameters 
especially  the  Flags,  is  available  m  the  problem  output  shown  in  Appendix  B.. 
The  input  is,  in  general,  standard  NAMELIST  format  (see  listing  of  sample 
data  in  Appendix  A)..  Following  the  NAMELIST  input  is  a  single  card  con¬ 
sisting  of  a  case  description  (any  80  alphanumeric  characters)  which  is 
printed  at  the  start  of  the  case.  To  signify  the  end  of  the  run  (no  more 
cases),  a  NAMELIST-format  input  is  required  setting  the  flag,  END,  equal 
to  1.  Mo  case-description  card  should  follow  here. 

A  sample  problem  follows  for  an  engine  typical  of  a  solid-propellant  inter¬ 
ceptor  sustainer  rocket.  Since  the  radiating  species  mole  fraction  (XM)  has 
been  set  equal  to  zero,  only  radiation  from  the  alumina  particles  will  show 
up.  The  specification  of  look  angle  equals  zero  implies  a  broadside  look. 

The  vertical  angles  equal  to  90  degrees  means  that  the  plume  axis  is  parallel 
to  the  Earth's  surface.  The  values  of  90  degrees  for  the  velocity  vector 
imply  that  the  free- stream  velocity  is  perpendicular  to  the  plume  axis.  The 
final  columns  of  printed  output  present  x,  y  coordinates  in  cm  followed  by 
the  intensity  of  w/cm^  sr;  many  zero  values  of  intensity  result  because  the 
particles  cannot  expand  through  the  total  plume..  The  output  from  the  sample 
case  follows:  (Note:  The  "error  summary"  in  the  output  refers  to  reading  in 
only  five  values  for  density  versus  time  instead  of  maximum  allowable  value 
of  50  —  no  calculation  errors  result). 

All  intermediate  calculations  le&ding  up  to  calculation  of  radiance  (emissivi- 
ties,  etc.,)  need  not  be  repeated  on  successive  cases,  if  not  changed..  These 
calculations  will  be  skipped  if  ISL  is  entered  as  -1  (solid  propellant)  or 
-2  (gas  only) 

Note:  If  ISL  =  +1  or  -1,  the  density  table  (PDE)  should  be  for  chamber 

density  and  not  exit  plane  density. 


Section  6 

RECOMMENDATIONS  FOR  FUTURE  WORK 

This  version  of  FLAME  forms  a  useful  framework  for  the  calculation  of  plume 
radiance  including  a  variety  of  mechanisms.  The  main  accomplishments  of 
this  contract  were  to  formulate  the  integration  procedure  and  geometric- 
orientation  choice  to  relieve  the  engineer  of  tedious  extrapolations  to  attempt 
to  relate  flight  data  to  calculations.:  Because  of  the  large  number  of  options 
included  in  FLAME,  a  comprehensive  checkout  was  not  possible  during  this 
contract  (over  50  runs  were  made).  A  streamlining  of  the  output  and  calcula¬ 
tion  procedures  would  be  possible  through  a  more  extensive  exercising  of  the 
FLAME  options.  A  variety  of  values  for  look  angle,  altitude,  times,  orienta¬ 
tions  engine  characteristics,  relative  velocity,  etc.  should  be  used  as 
FLA.vir,  input  to  see  if  further  refinements  are  needed. 

Currently  the  code  calculates  only  half  of  the  plume  in  a  single  run  which  is 
sufficient  for  axisymmetry,  such  as  for  particle  radiation  only.  The  code 
should  be  refined  to  allow  total  plume  extent  to  be  calculated. 

A  better  contour -plotting  routine  might  be  found  by  limiting  the  closeness  of 
the  isointensity  lines,  and  constraining  the  curve  fit  to  obtain  smoother  data. 
Automatic  labeling  of  the  contours  would  be  a  nice  added  feature. 

The  theoretical  models  used  in  FLAME  are  becoming  rapidly  antiquated.  A 
molecular  radiation  model  based  on  non- Boltzmann  distributions  should  be 
investigated  and  incorporated. 

Additional  extensive  flow-field  revision  would  be  useful  based  on  sophisticated 
multiphase  transient  plume  codes, 
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Appendix  A 

LISTING  OF  FLAME  CODE 


c  PROGRAM  EL  A  Ml  ,  VERSION  l 

PROGRAM  P?1 70(  1 MPUT, OUTPUT, F  ILMPL  ,TAPF5=  INPUT  ,  T  A  PE  6.  =OUT  PUT  , 

1  TAPF4R  =  FRMPL  ,1  APE16  ) 

DIMENSION  PX( 3),PY{  1  00) ,PL ( 3,100) ,FC (20 ) , IP2 (20  1 

0 1  MENS  ION  PDF  (  SO  ) ,  T  1  ()(  SO  )  ,  WL  E  (  ft  )  ,  F  PT  (  ?  0 , 1  0  )  ,  T  FT  (  2  0  )  ,TCOOL  ( 

120,10)  ,TC.1  (20)  ,BPASS(  20)  ,FTL(  20,10)  ,TMIN(  10  ) ,  DEL  T  (10)  ,ET  (20),  1  R(  20 
2 ) , JM(20) ,JT!20) ,TJ  (20) ,XP! P) ,W(  8) ,RP ( 10) ,TC00(20) ,TC 1X120, 10), CASE 
3<  10)  ,T  I  0P(  51  )  ,,  PDF  P  (  SI  )  ,T1  (  10)  ,  T2  (  10  )  ,WP  (20  )  ,CU(2  ) 

EOtll  VALENCE  (TC  IX.T10P)  ,  (TC.IXII  ,  4  )  ,  P  0  E  P  )  ,(T1,JT),(T2,JT(111) 

REAL  IA,MF,.)M,JT,L1P 
INTFGFR  FNO 

COMMON/ MCG/  JO, I  PL , PX , PL, PY,FC ,XL ,XR,YR ,YT , IP, JP2 
COM M. 0 N / M S T C F /  !  A  ,  Al.  ,  RO  ,CONF  ,RM  ,  X 1  ,Y1  ,  R  A 

COMMON/ MF/ 1 SL  «  NP, G*  RE  *  GAMMA ,ME , AN ,T I G,T I D, POE  »NR , XML  *OELT , TMIN, TCI 

1  ,  TC  OOL ,PI,RP,1E1 ,ETL,RPASS,TF,VP,P2,  X  1 2 , TR , ET , A , VI , V2 , V3 , HO *  WMOL 
2, XM, XP , W, GL 

OAT A  XP , W , R A ,P1 ,RD,WP/ -.9602 898565, -.7966664774, -.52 5532 4099, -.183 

14346425. .  1834346425. . 5255324099. . 7966664774 . . 9602 898 565 . . 101228536 

2  3, .2 2? 38 10344,. 31  37066458 ,. 3626837834 , . 362 68 37 834 , . 31 37 0664 5 8 , . 2 22 

33810344. .  1012285363.90. .3. 141 5927.57.29577951 31 . .001 .1 9*1. / 

OAT  A  NP,RP/7,46.S04E-5,40.E-5,3?.-31  5F-5,26.6E-5,?0.001E-5,13.3015E 
1 -5, 6. 9007 F-5/ ,8/9.49/ 

OAT  A  TR/ 100., 400., 800. , 1 2 00 . , 1 600 . , 2000 . ,2 400., 2 700. ,3000., 3300. , 

13600. ,  3900. ,4200. ,4500. , 5000. , 5 500 . , 6000 . , 7000 . ,9000. ,11000./ 

DATA  TF1 /99.  ,1 50. ,200. ,2 50., 300. ,350. ,400. ,45u. ,500. ,600. ,700. , 

1R00.,  1000. ,1200. ,1400., 1700., 2  000., 2 500., 32 00., 6001  ./, 

2  TF.OF, CP/2323. ,4240. ,3. 2/ 

C  *****  REGION  0,  INPUT  OATA  ***** 

NAMELIST  /IN/Tli,  ,RE,GAMMA,ME,L1P,TE,AL,H0,VR,PHIU,XIU,TC,PHIR,XIR, 
1 WLMIN,WLMAX,XM, I  NO, ISL, 1NT,LINES,NR,TID,PDF,END, I  PLOT 
GL  =  A LOG! 10 . ) 

IPT  =  0 
1  READ! 5, IN) 

IF  (FNO. E 0.1 )  STOP 
JP  =  1 
KP  =  1 
I  PL  =  0 

00  32  I  =  1,20 
32  IP? (  1  )  =  0 
IP  =  0 

READ (5,200)  CASE 

WRITE (6, 100)  CASE,TIG,RE,GAMMA,ME,L1P,TE,AL , HO , VR , PHI  0 , X I U  ,  TC  ,  PH  I  R 

1  ,X!R ,WLM1N,WLMAX,XM,IND, ISL , I  NT .LINES, NR, I  PLOT 
200  F0PMA1  ( 1 0 A 8 ) 

100  FORMAT ( IH1 , 3 X ,6 7HPR0GR AM  P2170  FLAME  ROCKET  EXHAUST  PLUME 

10PT1CAL  SIGNATURES  / / /POX ,  1  0A8//1 1H  INPUT  DATA//29H  TIME  AFTER  IGN 

2  IT  1 ON-T  !G(SEC  )F10.4,11X,2  5HN0ZZI.  E  EXIT  R  AOI US-RE  ( CM  )  FI  0.4 , 13X  ,  5HG 

3A MMAF 10.4/  X,19HEXIT  MACH  NUMBER-MF F 1 0 .4 , 2 OX ,2 5HN0ZZL E  LIP  ANGLE- 
4L I P ( OEG ) FI 0 .4 ,1 3X ,26HEFFECT I VE  FARTH  TEMP-TE (K ) F8.2/  X,18HL00K  AN 

5GLF-AL ( OEG ) FH . 2 .2 3X,22HENGINE  ALT ITUOE-HO ( FT IF9.0, l 7X ,23HROCKET  VE 
6L0CITY-VR(FPS)F9.2/  X,27HVERT.  POLAR  ANGL E-PH 1 U ( DEG ) F8 .2 , 14X , 30HV 
71 RT.  AZIMUTHAL  ANGL E-X IU ( OFG ) F 8.2 , 1  OX , 1 9HCHAMBER  TEMP . -TC ( K ) E8 . 2 / 

8  X  ,2  6HVEL  .  PULAP.  ANGLE-PHI  R  (  OEG  )  E8  .2 , 1 5X  ,  29HVEL  .  A7  IMUTHAL  ANGLE-X 
9IR(DEG)FR.2/16H  BAND  LIMITS  AREF6.2,3H  T0F6 .2 , 1 9X , 30HR ADI  AT  I NG  SPE 
ECUS  MASS  FRAC-XM  F8.5,//  23H  CONTROL  FLAGS  F 

AOLLOW — / / 8H  IND  =I5,25H  (1  FOR  DAY,  0  FOR  NIGHTI/8H  ISL  =15,79 

BH  (1,-1  FOR  PARTICLES;  0,-2  FOR  GAS  ONLY;  NEGATIVE  TO  SKIP  INITIAL 
G  CALCULATIONS)  /8H  INT  =I5,56H  (APPROXIMATE 

CNUMRER  OF  INTEGRATION  STEPS  THROUGH  LAYER ) /8H  LINES  =I5,56H  (APPRO 


o  o  o 


DXIMA1E  NO  ."HER  OF  L  1  NF  S  OF-  SIGHT  IN  ONE  DIMENSION  )/8H  NR  =IS,44H 
F  (NUMBER  OF  FNTRIFS  IN  FIRING  HISTORY  TABLFI/BH  l^LUT  =!S,42H  (0  F 
FOR  NO  PLOT,  1  FOR  ISOIMTFNSITY  PLOTS)) 

C  *4444  INITIALIZATION  OF  SO  40? 0  CAMFRA 
IFF  I  PLOT  .  ME  .<> )  GFJ  TO  30 
31  IR  =  (NR+?)/3 
K  =  l 

DO  2  1  =  1  ,  I R 
00  ?  J  =1,3 
L  =  I  -r  (J-l  )=?IR 
TIDP(K)  =  TID(L) 

PO  F  P  (  K  )  =  PDF  (I.  ) 

2  K  =  K  +1 

WRI1F  (6,1 01  )  (T I  OP  <  1  )  .PDFPtl)  ,1  =1,NR) 

101  F  OR  HA  T ( 7 2 HO  FIRING  HISTORY  TARlF  FOLLOWS,  EXIT-PLANF  DENSITY  (G/CM 
13)  VS  T  1  ME  (SEC)  //3( 4 X, 4HTI Mfc,5X,7rli)FMSITY, 5 X)/(3(F10. 4,611,3,4X1 
2  ) 

#4***  REGION  1,  INITIAL  CALCULATIONS  -  PARTICLES  44444 

CC  = ( WLMfl  X  -  WLMIN)/?. 

00  = ( WLMA  X  +  WLMIN  )/2. 

I F ( ISL.NF.l )  GO  TO  3 
DO  4  1  =  1,2  0 
T  =  TE1 ( 1  ) 

C  FIND  WAVELENGTH  WHERE  PLANCK  FUNCTION  IS  GREATEST.  TRANSFORM  INTEGRAL 
WW  =  2R98./T 
C  =  3.94WW 
0  =  4.14WW 

C  PFRFORM  GAUSSIAN  INTEGRATION  (EQUATIONS  4-1 , -2 , -4 , -5 ) 

BRASS !  I  )  =0. 

Rl X  =  0. 

00  7  K  =  1 ,NP 
EPT (  I  ,K  )  =  0, 

7  ETL ( I  ,K  )  =  0. 

00  5  J  =  1,8 
Wl  =  C4XP ( J  )  +  0 

wwl  =  cc*xp(.i  i  +  on 

Rl  =  3.742E4/WL44S/ (EXP( 14390./WL/T 1-1 . I 
R2  =  3.742E4/wWt4*5/(EXP(14390./WWl/T)-l.) 

Rl X  =  Rl X  +  W ( J ) 4Ri 

BRASS ( I )  =  BPASS(I)  +  R24W(J) 

C  SUM  OVER  PARTICLES 
00  5  K  =  1 ,NP 
EMI  =  EMIT (RP(K),WL) 

EM2  =  EMIT (RP(K ) ,WWL I 

EPT( I ,K)  =  E  P  T ( I , K  )  +  WIJJ4R14EM1 

5  ETL (  I  ,K I  =  EIL(I.K)  +  WUI4R24RM2 
00  6  K  =  l, HP 

EPT (  I  ,K )  =  FPT ( I ,K  I/R1X 

6  ETL ( I ,K I  =  ETL ( I , K ) /B PASS (  I  ! 

4  RPASS(l)  =  BPASS( I I4CC/3. 14159 
WR I TE ( 6 , 1 03 )  ( RP( I ) » 1=1 »NP) 

103  FORMAT ( 46H0  INITIAL  CALCULATIONS  FOR  PARTICLE  PROPERT I ES //35X, 16HT 
10TAL  FMISS1V1TY//16X,2 7HPART ICLE  RADIUS  (CM  X  10E4I/13H  TEMP ( K ) 

2  .4P10F10.3//) 

00  20  I  =  1 ,20 

20  WR 1 T F ( A , 104 )  TET ( I  I , ( E PT (  I ,  J I , J =1 , NP ) 

104  FORMA  T ( 1H  , F9 . 1 , 3 X , 1 0F1 0 . 5 ) 
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WR  IT  F (6, 107 )  WLMIN.WLMAX , (RP( 1 ) , 1=1 ,NP ) 

107  FORMA i (////] 3H0  RAND  DATA, ( ,  c  5 . 1  ,  4H  TO  , F 5 . 1  , 1  H ) / / 32 X , 1 RHFM I  SSI  V I  T 
W  IN  RAN!)//]  2  X,  42HPI.  ANCK  FUNC.  PARTICLE  RADIUS  (CM  X  10E41/25H 

?  TFMP(K)  (W/CM2/SR)  .4P10F9.3//) 

OH  ? 2  I  -  1  ,?0 

?2  WRITF  (f>,108)  TFT<I),RPASS(I>,(FTL(I,J),J  =  1,NP) 

108  FORMAT  (1H  ,  F9.1  , 1'l  3 .4,2X  ,  1  0F9.5  ) 

C  COOL  IMF,  CURVES  FOR  SOLID  PART  ICLFS  (SECTION  4.1.?, EOS.  4-3  AND  4-8) 

DT  =  ( TF-100 . )/ 1 9. 

TCOI)  (  1  )  =  TF 
DO  8  K  =  1 , NP 
TCIV( 1 ,K)  =  0. 

Tl(K)  =1  ,/TF**4*RP(K)  / T  ARC  E ( TF.TET ,EPT(1,K),20,0,1 ) 

8  DEI. T(  K  )  =  T1  !K  J-0F/..7016E-11 
T  =  TF 

DO  9  I  =2,?o 
T  =  T  -  DT 
tcuo ( i )  =  r 

DO  9  K  =  1  ,  N  P 

T2(K)  =  1 ./T*~4*RP(K )  /TARLFIT, TFT, EPT(1»K), 20,0,1  ) 

TCIX(1,K)  =  TC  I  X  (  I -1  ,K  )  +  (  T?  t  K  )  +  T 1  (  K  )  )  /2  •  *DT 

9  T  1  ( K  )  =  T2(K  ) 

DO  10  K  =  1 ,NP 

10  TCI X( 20, K  )  =  TC I  X ( 20 , 1  ) 

C  RFARRANGF  TO  CHANCE  independent  variarle  to  time 
OT  =  TC I  X ( 20 , 1  J/722700. 

TIME  =  -DT 
DO  11  I=i ,?o 

TI^E  =  TIME  +  DT*FLOAT ( 1**4 ) 

TCI ( 1 )  =  TIME/1. 7016E— 11 #CP 
DO  11  K  =  1 , NP 

11  TCOOL(I,K)  =  TARLE(TIME,TCIX(l,K),TCOO,2O,0,l) 

TC  I  (  ?0  )  -  1 . E5 

WRITE (6il05)  (RP(  I  ),  I=lfNP) 

105  FORMAT (/ //2 b X ,34HPART ICLE  TEMPERATURE  HISTORIES  ( K ) / / 1 6X , 2 7HP AR ( 1 C 

1LE  RADIUS  (CM  X  10F4J/13H  T I  ME ( SEC )  , 4P1 0  FI  0 . 3// ) 

00  21  I  =  1,20 

21  WR  I  TF  (  6 , 1  Of> )  TCI  (  I  ),  (TCOOLC  I,J  ),J  =  1,NP) 

106  FORMAT  ( 1  FI  ,  F 1  0 . 3  ,2  X  , )  OF  1 0 . 1  ) 

3  IF(  ISL.LE.O)  GO  TO  1'. 

TO  =  400. 

DO  12  K  =  1 , NP 
TC  00l ( 2  0 ,K )  =  70. 

C  MEWTON-RAPHSON  SOLUTION  TO  FO.  4-7  (AND  EO.  4-6) 

C  =  5.964E10*TARLF(6000.,TET,EPT( 1 ,K ) ,20,0,1 )#FLOAT( IND)  + 

1TARLF  (  TE,  TET  ,EPT  (  1  ,K  )  ,20,0,1  )!*TE**4/2. 

DO  13  1  =  1,100 

TN  =  (  2  .  $  TO  +  C  /T04-«3/TABLE(T0,TET,EPT  ( 1  ,K  1,20,0,1  )  1/3. 

IE  ( A8S( TN-T01.LT..1 )  GO  TO  14 

13  TO  =  IN 

WRITE  (6,102 )  TN , TO ,K 

102  FORMATf S8H0ERR0R,  SOLUTION  OF  FQ.(4-7)  DOES  NOT  CONVERGE.  NEW  VALU 
1E=  G8.2,  12H,  OLD  VALUE=  G8.2 , I  5 , 1 RH-TH  PARTICLE  GROUP) 

14  TMIN(K)  =  TN 

12  CONTINUE 

WRITE(6,109)  ( RP  (  I ) , I  =  1,NP) 

109  FORMAT ( / // /2  7H  PARTICLE  RADIUS  (CMX10E4)  ,4P10F10.3) 

WRI TF ( 6, 1 1 0 )  ( DEL  T ( I ) , 1  =  1,NP) 

110  FORMAT ( 2  7H  SOLIDIFICATION  TIME  (SFC)  .10F10.3) 


u 
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URI TF (6,1 1 1 >  (  T  M I  M  (  I  )  ,  I  =  1  ,NP  ) 

111  F 11 R  M  A 1  ( 2  7  H  F  01 1 1  L  1  BR  HIM  T  fcMPERA  CUR  E  (  K  )  .10F10.2) 

GAS 

RUN,  SKIP  RANO-FNfcRGY  CALC. 

1  6  IH  ISL.l'  .0)  GO  !  0  16 
DETERMINE-  FI  (FO.  A  — 1 1  ) 

Ul  =  0.5FA/R/WLMAX 

o?  =  o . sfa/r/wln in 
on  ir  l  =  i ,?n 
t  =  T  k  (  n 

S  TG  -■  1  .43*7-6  /T 
XMIN  =  01  *01  ~-S  I  G 

x  m  a  x  =  r 

FT  (  !  )  =  9.928E-4*R/SIG$(01-XFXf>(-XMJN)  -  0?  *6  X  P  (  -XM  AX 

11  h  . MR62 3X30RT( S IGI* ( FRF ( 02*SGR1 ( S 1G  )  )  -  F RF ( 01 -SORT  ( S I G  1  )  )  ) 

OFT F KM  INF  BA  NO  RADIANCE  (ED.  4-15) 

JM( 1 )  =  0. 

DC)  1  7  J  i  ,  B 

wwl  =.  cc-xp(j)  +  on 

17  JM( I  )  =  JM(I>  +  W( J ) SAHCOf WWL ,T ) /WWL**5/ ( EXP (1  A 390. /WWL /T ) -1 . ) 

18  JM(I)  =  JM(  I  )  aCC*3.742F4/P1 
WRIIFIA.IP)  (TR(1),ET(1),JM(1),I  =  1,20) 

112  FORMA f ( ////A3H0R0TAT IONAL  RADIATION  FROM  WATER  MOL FCUL ES . /97H 
’.AVAILABLE  RAC'D  FNERGV  =  ET  (10E-20  W  SFC/MOLFC),  BAND  RADIATIVE  RA 
2TE  =  J  ( 1  OF -2,0  W/  SR  MOL  EC  )  //?4X , 4HT (K ) ,8X , 2HET , 12 X , 1 H J// ( 22 X , F7 .  1 
3,2  K,?HUI  1 

CALC0LA110N  OF  RA0IA1IVE  COOLDOWN  ( fc 0 .  4-16)  (J  VS  TIME) 

ET  X  =  e T ( 1  ) 

TO  =  0. 

DEI  =  CT (201/190. 

.1  T  (  1  )  =  j  M  (  1  ) 

TJ ( 1  )  =  0. 

DO  19  I  =  2,20 
ET  X  c  ET  >•  +DET 

J  T ( I )  =  TARLE< ETX- CT , JM, 20 , 0, 1  ) 

TN  --  .0/9r'7R/JT  1  I  ) 

T  J  (  I  )  =  TJC-I)  -(TO  +  TN1/2.-DET 

19  TO  =  TN 

FIT  DATA  WITH  FUNCTION  LINEAR  IN  TIME  (EO.  4-17) 

CALL  LSOPOL  (TJ,.IT,WP,T1DP»20»ETX,2  ,CU) 

A  =  C(J  (2  ) 

DO  23  1=  1  ,20 

TCOOI I )  =  CD( 1 )  +  CU(2)*1J( I) 

T1DPU)  =  TABLE  (  JT  (  I  )  ,JM,ET,20,0,1  ) 

23  PDFP(I)  =  TABLE!. I T( I  ),JM,TR,?0,0,-1) 

WRITE!  6, 113)  (T.l '  I)  ,JT(  1  )  ,TCOO(  I)  ,TIDP(  1  )  ,PDEP(  I )  ,  I  =1,20) 

113  FORMAT! ////28H  RADIATIVE  DECAY  IN  RAND  //9X,  9HT I  ME ( SEC  ) , 6X , 1 H 

1J,BX,6HJ(FIT) ,7X,?HET,8X,4HT(K)//(6X,F11.5,2X,4E11.3) ) 

SOLUTION  FOR  LIMITING  VELOCITY  (EOS.  4-9  AND  4-10) 

16  VP  =  2 .-SORT ( TC /NMOL / ( 2 • -2  ./GAMMA  )*R.3144E7) 


region  2,  INITIAL  CALCULATIONS  - 
IF  BAND  LIMITS  UNCHANGED  FROM  PREVIOUS 


PREPARATION  FOR  PLUME  INTEGRATION 

XI  =  0. 

Y 1  =  0. 

C  CALCULATION  OF  MAXIMUM  PLUME  CONE  HALF-ANGLE  AMD  MAXIMUM  PLUME  RADIUS 
GG  =  SORT ( (GAMMA+1 . )/ (GAMMA-l . ) ) 
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.  I  -- .  ■  »>  "..'I'-MiMWie 
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(>M  =  S0RT(MF**2  -  1.) 

conf  =  up  +  ra*igg-ii  -  rd*(gg*atan(gm/ggi  -atan(gki> 

RM  =  VP*TIG 

DELS  =  (RM  +  RM*COS<  AL/RD  )  )*.05/  (  1  .1  **(  1NT/2  )  -  1.) 

IA  =  n. 

PI  =  PHIR/RI) 

P2  =  PHIU/RI) 

XU  =  XIR/RO 
XI?  =  XIU/RD 
IF( AL.6T.45. )  1A  =  1. 

DEL  X  =  DFLS*COS  (  <RA*IA  -  AU/RD) 

DEI.  Y  =  (RM  +  RM*SIN(AL/RD>  1  * . 1 / ( 1 . 1 **L I  NFS  -  1.) 

CD  =  COS  (  (  RA*1A-AL  )/RO) 

OELXX  =  DEIY/CO 
DEL  =  DELY 
JO  =  10000 

C  SELECTION  OF  1S0INTETS1TY  LEVELS  TO  BE  PLOTTED 
FC( 1 )  =  1  - E-l 1 
XI 0  =  SORT (IP.  ) 

DO  29  I  =  1,19 
?  9  FC(  1+1  1  =  FC ( 1 ) ’ X 1 0 

*****  REGI ON  3.  SELECTION  OF  INITIAL  VALUE  OF  XI  ***** 

VV  =  VR*30.4801 
VI  *  V V * S 1 N ( P 1  )*COS<  XII ) 

V?  =  VV*S I N ( PI  1 *S I N ( X  I  1  ) 

V  3  =  VV*C0S(P1  ) 

IN  (JUT  =  ? 

I F ( 1A.EC.1.)  GO  TO  36 
1 F ( XM.EO.O. )  GO  TO  27 
I F ( CONE .GT . 90 . )  GO  TO  25 
1  F  (  90  .  -  CONE.GF.AL)  GO  TO  27 

34  XI  =  RM*COS ( ( C  ONE  +  A L ) / RD 1 /C OS ( AL/RD ) 

GO  TO  27 

25  I F ( 90  .  -  CONE.GE.AL  -  90.)  GO  TO  34 
XI  =  -RM/COS( AL/RD) 

GO  TO  27 

36  I F ( XM.Nfi .0. )  GO  TO  33 

XI  =  -RM*COS ( ( AL  -  33.52 )/RD)/SIN( AL/RD) 

GO  TO  27 

33  IF(CCNF.LT. 90.. AND. CO NE.LT.AU  GO  TO  35 
XI  =  -RM/S  I N  (  AL  / Pi) ) 

GO  TO  27 

35  XI  =  -RM *C OS i ( AL  -  C ON E)/RD)/SIN (AL/RD) 

27  XPLOT  =  0. 

PX( 1 )  =  XPLOT 

DEL  XI  =  (OELXX  +  ABS(X1)/10.)/1.1 
R  =  1 ./I  .1 

PLOT  LIMITS  FOR  GRID 
XL  =  0. 

XR  =  RM  -  XI  *COS  (  (  R  A  *  I A  -  AU/RD) 

IF( IA.EO.l  • .AND .XM.EO.O. )  XR  =  2 . *RM*S 1 N ( 33. 52 /RD ) 

1  *S I N( AL/RO ) 

YB  =  0. 

YT  =  XR 

*****  REGION  4.  CALCULATE  INTEGRATION  BOUNDARIES  ***** 
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2 4  CALL  STAR(°M1N,PMAX, IN0U1  I  ’ 

OFL  V  =  DF.LY*1.1 
I  F  (  1NOUT.FO.O)  GO  TO  2  6 
Y1  =  0. 

nr-LY  =  of  l 

JO  =  M 1  NO (  J  0  ,  J  P  I 

KP  =  KP  +  1 

IF(KP.EQ.4.AN0.  IPLOT.NF.O)  CALL  CONTUR(KP) 

Jp  =  0 

IF  ( XI .GF.O. )  R  =  1.; 

0ELX1  =  OFL  XI  *R 
XI  =  XI  +  OFL  XI 
XPLOI  =  XPLOT  +  DELX1*CD 
C  STORE  VALUES  (PX.PL.PY)  for  plotting 
PX(KP )=XPLOT 

CALL  STAB(PMIU,PhAX, INOUT  ) 

I  F  (  IUOUT.fcO.O  >  r,0  TO  26 
GO  TO  l 

*****  REGION  5,  INITIAL  EVALUATION  OF  INTEGRAND  ***** 
26  CALL  FUN ( SOL 0  » PM  IN ) 

*****  REGION  6,  INTEGRATE  THROUGH  PLUME  ***** 

CALL  S IMP (SOLO, SI  NT, PM  AX, PM  IN, DEL  X  ) 

SINT  =  SINT*. 6666667 
WR IT  F ( 6 , 1 1 6  I  XPLOT ,Y1 , SINT 
116  FORMAT(4X,2F10.2,E11 .31 
JP  =  JP  +  1 
PL(KP,JP)=  SINT 
PY(JP)  =  Y1 
Y1  =  Y1  +  DELY 
GO  TO  24 

30  IF(IPT.EO.O)  CALL  C AMRAVt I  PLOT ) 

IPT  =  1 
GO  10  31 
END 
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SURKOtn  INF  S IMP <SOLO,S I  NT , PMAX , PM  I N , OFL ) 

sint  =  smn 

OELX  =  (  ARS(PMIN)  /10.  +  OELl/1.1 
A  =  1  ./I  .1 
X  =  PM  IN 
1  SUM  =  0. 
f.  =  4. 

OELX  =  ORLXSA 
A  X  =  X  +  OELX 

I R( X.GT.PMAX)  GO  TO  ? 

CALL  FHNtSS.X) 

SS  -  SS*OELX 
SUM  =  SUM  +  C* SS 
IFIC.EO.a. )  GO  10  3 
I  E ( X.GT.O. )  GO  TO  b 

5  SOLO  =  SS 

SINT  =  SINT  +  SUM 
GO  TO  1 
3  C  =  ?. 

GO  TO  A 

?  SINT  =  SINT  -  SOLD 
RETURN 

6  A  =  1 .1 
GO  TO  5 
END 
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SUBROUTINE  STAR( PMIN, PMAX, INOUT ) 

DIMENSION  X  C  2 ) 

COMMnN/STC/Xt  XX 

COMDON/MSTCF/  I A  ,  AL  ,RD,CONE ,RM,X1  ,Y1 ,RA 
REAI.  !  A 

1 F ( INOUT.NE.2 )  GO  TO  1 

AO  =  ABS  <TAN( (RA-IA  -  AL1/RD )) 

Al  =  AO*AO 

A?  =  ( TAN< < RA*( I A-l . CONEJ/RO) )**2 

1  INOOT  =  0 

R1  =  XI *AO/ ( 1  .  +  Al  ) 

R  =  Y1SS2 

Cl  =  B  +  XI --2  -RM*RM 

1 F ( 1 A.EO.O. )R  =  -B~A2 

C2  =  XI Y  XI  +  R 

RAD  -  B1  *R1  -0/(1.  ■*  Al  ) 

J  I  =  1 

C  CHECK  FOR  L INE-OF-S IGHT  BEYOND  R=RP. 

1 F( RAD.LT. 0. )  GO  TO  5 
RAD  =  SORT (RAO) 

XX  =  -B 1  +  RAD 
CALL  CHECK (II) 

I F( II . GT . 2 )  GO  TO  4 
XX  =  -HI  -  RAD 
CALL  CHECK ( I  I  ! 

I  F  (  II.GT.2)  GO  TO  4 

3  RO  =  Al  -  A? 

R2  =  X1*A0/BD 

C  CHECK  FOR  LINE-OF-SIGHT  BEING  PARALLEL  TO  PLUME  BOUNDARY. 
IF( ARS(BD) .LT.l ,E-B)  GO  TO  2 
RAO  =  B2*B2  -  C2/BD 
IF  (RAD.LT.O.)  GO  TO  5 
RAD  =  SORT ( RAD ) 

XX  =  -R2  +  RAD 
CALL  CHECK ( I  I  ) 

1F( 1  I.GT.2 )  GO  TO  4 
XX  =  -B?  -  RAD 
7  CALL  CHECK (II) 

I F ( 1 1 .  G  T  .  2 )  GO  TO  4 

5  INOUT  =  1 
RETURN 

2  IF(AL.EQ.O..OR.AL.EQ.RA)  GO  TO  5 
XX  =  C2/R2 

GO  TO  7 

4  IF  ( X( 1 ) .GT.  X (2 ) )  GO  TO  6 
PMIN  =  X(l) 

PMA  X  =  X  (  2  ) 

RETURN 

6  PMIN  =  X(2) 

PMA  X  =  X ( 1  > 

RETURN 

END 
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SURROUT  IMF  CHECK  (1) 

D1 MENS  I  DM  X(2) 

COMMON/STC/ X,  XX 

COMMON/ MS  TCP  /  1  A,  AL.Rl),CONE,RM,  XI  ,Y1  ,  RA 
REAL  I A 

C.  CHECK  TO  SEt  IP  XX  REPRESEN1S  X  OR  Z. 

1F(  I  A. E  O.O.)  (,0  TO  1 
17  -  XX 

RR  =  SORT  (  Y1  =>YT  +  t  Z  Z  *T  AN  I  (  R  A-AL  ) /RO )  +  Xl)**2  +ZZ*ZZ> 

2  IP  (RR.EQ.O.)  RETURN 
THETA  =  RD«ACOS l 7Z/RR) 

C  CHFCK  TO  SEE  IF  INTERSECTION  IS  ON  CONICAL  PLUMF  ROUNDARY. 

IF ( ARK  I ( THETA -CONE  1/CONE  )  .LT  *.001  )  CO  TO  3 
C  INTERSECTION  IS  on  SPHFR1CAL  CAP  OF  PUJME. 

IF  I  THETA .GT .CONE  I  RETURN 

C  00  THF  FOLLOWING  IF  THIS  POINT  IS  INDEED  ONE  OF  THE  INTEGRATION  LIMITS 
4  XII)  =  XX 
I  =  1+1 
RETURN 

1  77  =  XX*TANt AL/RD)  +  XI 

RR  =  SORT ( Y1 *Y1  +  XX*XX  +  ZZ*ZZ) 

GO  TO  2 

C  INTERSECTION  IS  ON  CONICAL  ROUNDARY 

3  I P( RR/RM.LT.l .001  I  GO  TO  4 
RETURN 

END 
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SURROUI INF  FUN (0, P) 

DIMENSION  PDF  (50)  ,T  1 0  (  50)  ,  TFT  (20)  ,  TCOOl.  (  20 , 10  )  , TCI ( 20 ) , B P ASS  I  20 ) , 
IF  Tl  (20,10  )  ,TM1  N<  10)  ,DELT(  10)  ,  FT(  20)  ,TR  <20  )  ,RP(10  ),N<  10),  ARUO)  ,B(  1 
20) ,0  HO) ,TH(10) ,D( 10) ,E< 10),F(10) ,H1 (20) ,TMT(20> , CAT (20) ,XP(8> 

3 ,  W  (  n  ) 

RF.A1  ME,  IA  ,J  ,.IG,N 

COMMON/ MF/ISl , NP,R,PE,GAMMA,MF ,AN,T 1 G , T I  0 , PDF , NR , XML , DFL T , TM I N , TC I 

1  ,TCOUl. ,  J  1  ,RP,1F1  ,FTI.  ,RPASS,TF.VP,PH1U,XIU,TR,ET  ,A,V1  .  V2  , '<3 , HO,  WMOL 
,  XM ,  x  P ,  W ,  Gl 

COCMON/MSTCF/  IA  ,  Al_  *  RD.CONF  ,RM ,  X  l  ,Y1 ,RA 

DATA  WM0L.XML.0P/17.  ,  .35  , A  .  0045 / , S  ,C F / .  7F-1 5 , 3 . 531 445E -5 / 

C  C.URVF  FIT  C0FFF1C1FMT  V  AND  RHO  CORRELATIONS  FOR  TWO-PHASfc  FLOW 

DATA  AB/2 . 500  7 F 5, 2. 676405,2 .7520F5.2 .H248E5,2 . 89 32 E 5, 2 .97  OIF  5, 3. 07 
175F5  /  ,0/-2  . 9 7 A 5 F 5  ,  — 2  .9 337E4  , -2 . 935  5E 4  , -3 . 0367 E4 , -3 . 2  891  E4  , -3 . 8094F. 

2  A  ,- A.  5592  F'»/,C/.  02092  ,  .0205A,  .  02  035  ,  .02  01  2  ,  .0 1  97  6  ,  .0 1  930,  .  0  187  A  /  ,  N 
3/-. 5 A? 5,-1 .0308,-1 .31 75,-1 .4886,-1 . 5A 1 9 , -1 . 349 A , - . 9048/ , 0 / . 5 395 , . 8 
a  87, 1 . 3339,1 .0308, .5395, .4053, .2201/, E/2. 04, 2. 05, 2. 08, 2. 12, 2. 17, 2. 2 
54,2 .33/, F/-2 .05, 5. 30, 5. 95, 7. 30, 5. 75, 5. 55, 2.84/, TH/1  3.56, 14.7  5, 15.8 
5 1,1 9. 76, 23. 53, ? R. 22, 33. 52/ 

C  H( FT ) , T ,  AND  LOG ( N/ F T3 )  FROM  IJ.S. STANDARD  ATMOSPHERF  (1962  ) 

DATA  HT ,TMT , CAT/3. F5, 3. 5F5,4.E5,4.5E5,5.F5,5.5F5,6.F5,6.5E5,7.E5,8 
1.E5,9.F5,1.F6,1.1F6,1.3F6,1.5E6,1.7E6,1.9E5,2.1F6,2.3E6,3.0E6,184. 
2  9  5,241  .74, 385. 3a, 663.20. 924. 15 ,1084. 09, 1 1 70.a7 , 1  2 30 , 3, 1  2 75 .2  1 ,  1346 

3.67. 1 396. 1 4. 1  A3  8. 68, 1 4 55. 53, 148 5. 69, 14R9. 35, 1495.2 9, 1500. 57, 1506.3 
4  3,2*1  507. 58,18.1 >,1 6.98, 1 6.05, 15.41  ,1 5.01 ,14.75, 14.5  5,14.36, 14.20, 

513.90.13.64.1  3.40,13.17,12  .76,12  .41 ,12.08,11 .76,1  1 .49,11  .20,  10.3/ 

X  =  P 

JG  =  0. 

SI  =  0. 

G2  =  GAMMA  -  1 . 

Z  =  XI  +  A8S(TAM( (RA*1A-Al)/RD))*X 
I F (  IA.E0.1 . )  GO  TO  6 
7  R  =  SORT! X* X  +  Y1$Y1  +  Z*Z) 

AN  =  ACOS ( Z / R )  *RD 
I F (  1  A RS ( I SL ) .NE.l  )  GO  TO  1 
Gl  =  (Gamma  +  1 .  )/2  . 

C  CALCULATE  THROAT  RADIUS  FROM  EXIT  RADIUS  AND  EXPANSION  RATIO 
RT  =  RE/S0RT(G1**(-G1/G2)*(1.+G2*ME*ME/2.)**G1/G2/MF) 

C.T=  COS  (  AN/RD  i 
DO  2  I  =  1 ,NP 

C  CALCULATE  PARTICLF  TIME  OF  FLIGHT,  ED.  4-25 
I F ( TH( 1 ) .LT.AN)  GO  TO  2 

C  ARE  WE  BEYOND  ThF.  LIMITING  STREAMLINE  EOR  THIS  SI7E  PARTICLE 
T 1  =  R/CT**N(  I  )/(AR( I )+  R( 1  )*FXP(-C( I )*R/RT ) ) 

I  F ( T I . G T . T I G )  GO  TO  2 

PDO  =  TARLFH IG-TI ,TID,PDE,NR,0,1  )*XML 
C  PARTICLE  CONCENTRATION 

PD  =  P00$0<  I)*(  RT/R  l«E(I)*CT«F(  I  ) /4. 1 8879/DP/RP  (  1  1**3 
C  CHECK  FDR  INCOMPLETE  SOLIDIFICATION  OF  PARTICLE 
I  F ( T  I  .LT  .DELT ( 1 ) )  GO  TO  4 

T  =  A  MAXI (TMIN( I ) ,TARL£(TI  -DFLT ( I ) , TC I , TCOOL ( 1 , I ) ,2 0, 0, 1 ) ) 

3  RO  =  TARLEtT, TFT, BRASS  ,20,0,1) 

C  PARTICLF  RADIANCE  ACCUMULATION,  EO.  4-27 

SI  =  SI  +  PO*PI*RP( I)**2  *rABLE(T,TET,ETL(l,I) ,20,0,-l)*B0 
2  CONTINUE 
Gtl  TO  1 

4  T  =  TF 
GO  TO  3 

C  GAS-PHASE  CALCULATIONS. 
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1  T 1 =  R /VP 

I F ( TI .GT.TlG.DR.AN.GE.CONF)  GO  TO  5 

GOO  =  TABLF<T[<;-T!,TID,PUE,NR,0,1  U6.0254E23/WMOL 

XK  =  OAMKA*G?  XMF*ME 

C  GAS  DENS  I T  Y  ,  FO  4-29  (MODIFIED  FOR  TWO-PHASF  FLOW) 

GO  =  GOO=!:X</?.*(Rfi/R,=x*2-XM4(COS(RA/CONE*AN/RD)  )*$(XK-2.  ) 

IF (GO  .LT.l.IGO  TO  5 
vx  =  X-VP/R 
VZ  =  7  *VP/ R 

C  RELATIVE  VELOCITY,  F.O  3-16 

VREL  =  SORT!  ( VX+V1  )**2  +  (  V7.  +V3  )  **2  +  (Y1*VP/R  +V2>**2) 

C  FIND  ALTITUDE  IN  PLUME  EO.  4-40 

H  =  HO  MS1N(PHIU)*(X*C0S(XIU  )+  Y1*S1N(XIU)I  +Z*COS( PH1U) ) /30.48 
H  =  AMAX1 (H,3.E5) 

H  =  AM [ N 1  ( H , 3 . E  6  ) 

C  FIND  EXCITATION  TEMPERATURE ,  FO.  4-35 
TA  =  TARLE(H,HT ,TMT,20,0,1  ) 

TM  =  AMAX1(TA,VREL*(9.61194E-9*VREL  +  1  .07636E-3  )  ) 

C  ATMOSPHERIC  CONCENTRATION 

CA  =( EXP( TABLE (H,HT, CAT, 20, 0,-1 )#GL ) I #CF 
C  CALCULATION  OF  ATTENUATION  FACTOR,  EO.  2-7 
HO  =  (H  -  H01/2. 

AT  =  0. 

00  B  I  =  1,8 

HS  =  HO  +  HD* (XP (I )  +  1  .  ) 

8  AT  =  AT  +  W( I l*EXP<TABLE<HS,HT,CAT,20,0,l)*GL ) 

AT  =  AT*R/2.“CP*S 

C  RADIATIVE  RELAXATION  TIME,  EO.  4-18  AND  4-36 
EM  =  TABLE(TM,TR,FT,20,0,1  I 
TAU  =  SORT ( EM/2. /PI/A ) 

C  COMPUTATION  OF  EMISSION,  EO.  4-38  AND  4-39 
TX  =  1.  -  TI/TAU 
J  =  A*TAU/2. 

I F ( TX.GT.O. I  J  =  J  +  J*TX 
C  RADIANCE  CALCULATION,  EO.  4-37 

JG  =  VREL*CA*S*TAU*(1.-EXPI-TI/TAU)  )*GD*J*1.E-20*EXP(-AT) 

5  0  =  SI  +  JG 
RETURN 

6  TX  =  2 
7  =  X 

X  =  TX 
GO  TO  7 
END 
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FUNC  1  ION  ARCO  (WL  ,TE  ) 

C  COMPUTES  ABSORPTION  COEFFICIENT  PER  MOLECULE  OF  WATER  IN  THE 
C  ROTATIONAL  SPECTRUM,  SCALED  BY  10E20 
1 F (  n .FO.TF)  on  TO  1 
T 1  =  T  F 
T  =  SORT ( T 1 ) 

WO  =  10.6-T 

YO  =  279.98/IT  SO. 031)  -  3.708 
A  =  1  .0303F-A  -  105.5R7E-4/IT-2.89M 

1  W  =  ] 0000 ./ WL 

ARCO  :  E  XP  (  7 .30?  5R8*(  YO  +  A*(W  -  WO )  /SORT  (  W )  )  )  *  1 . 3803E  -2  #T  E 

RETURN 

END 


FUNCTION  EMIT  (R  »W) 

C  C0MPU1FS  THE  FMlSSiVlTY  OF  ALUMINA  PARTICLES  OF  RADIUS  R  (CM)  A I  A 
C,  WAVFLFNGTH  W  (MICRONS)  KASFl)  ON  CORRELATION  OF  MIE  CALCULATIONS. 

I F  ( A.FQ.O. I  GO  TO  1 
IF  ( R.EO.RO)  GO  TO  2 

1  A  =  (11.32  -  2220. *R )  *R  +  .006R9 

RO  =  R 

2  EMIT  =  A*(W  +  2.S) 

RETURN 

El\!0 
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SUBROUTINE  CONTUR ( K  P  1 

n  I  MENS?  (IN  TT(9,XP(21  ),YP(2l),PX(3),PY(100),PL(3,100),FC(20),CA(6) 
1  ,  i  pr  ( ? n ) 

DAI  A  TT/9*10H  / 

COMMON/MCG/  ,.IO,  1  PI.  ,PX,PL  ,PY,FC,XL*XR,Yft,  YT,  IP,  1  P2 
JO  =11.11)  ^  1  )/2)*2  -  1 
! F  l  I  RL.FO.  0)  CALL  OR  1 D 
J1  =  JO  -  2 
FD  =  0. 

OX  =( PX(3)  -  PXIl I )/20. 

00  1  l.  =  1  ,  .11  ,  2 

C.  CURVE  FIT  TO  3  X  3  MATRIX  OF  PLUME  POINTS 
CALL  C  F  (  P  X  ,  P  Y  ( l_  )  ,PL()  »L  )  ,CA  ) 

C  HFRF  TO  STATEMENT  7  —  FINOS  WHAT  1  SO  I  NT  EN$.  1 T  Y  LINES  ARE  IN  MATRIX 
RIG  =  0. 

SMALL  =  1  .  F20 
L2  =  L+7 
00  2  I  =  1,3 
DO  2  J  =  L »L2 
A  =  RL (  I , J  ) 

I F ( A  .FO.O. )  GO  TO  2 
RIG  =  A  MAX] (R I G,A ) 

SMALI  =  AM  INI (SMALL  ,  A  ) 

2  CONTINUE 
VS  =  P  Y  (  L  ) 

YL  =  PY(L2) 

00  9  I  =  1  ,20 

IF  (PC( 1 >.LT. SMALL. OR. FC< 1 ).GT. BIG)  GO  TO  9 
IF  ( FO.GT.FCt 1  I)  GO  TO  7 
FO  =  FC(1 ) 

II)  =  I 
N  =  L 

C  HERE  TO  STATEMENT  3  —  SOLVES  FOR  COORDINATES  OF  ISOINTENSITY  LINES 
7  K  =  0 

X  =  PX<  1  ) 

00  3  .1  =  1,20 
R  =  CA<4)  +  CA(5)*X 
A  =  C  A  (  3  ) 

C  =  r,Alll*X«?  +  CA(2)*X  -  FCII)  +  C  A  (  6  ) 

0  =  R  =x=v2  -  4.*A*C 
IF  (O.LT.O.)  GO  TO  4 
0  =  SORT ( D ) 

Y  =  ( -R  +  0I/2./A 

IF(Y.GT.YS.AND.Y.LT.YL)  GO  TO  5 

Y  =  ( -R  -  01/2. /A 

I F ( Y.GT.YS.ANO.Y.LT. YL  )  GO  TO  5 
4  X  =  X  +  OX 

3  CONTINUE 

IF  (K.LE.l)  GO  TO  9 

CALL  AICRT3(0,0,XP,YP,K,1,1,1,5B,TT,TT,TT,2,1,1G.,16.,2,XL,XR,2, 

1 YR ,YT) 

9  CONTINUE 
1  CONTINUE 

00  6  M  =  l,  IP 

6  1 F (  IP2(M).F0.10)  GO  TO  10 
IP  =  IP  +  1 
IP2 (  IP)  =  10 

WRITE  (6,100)  EO  ,PX(1),PY(N) 

100  FORMAT ( 22H  ISOlNTENSITY  LINE  FOR  E9.? , 18HW/CM2/SR  IS  NEAR  (E9.2,1H 
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1 ,E9.? , 1 H  )  ) 

C  INI  T  Iftl.I  7ATIC1N  FOR  NFXT  GROUP  OF  3X3  MATRICES 
10  PX(1)  =  P  X ( 3 ) 
no  8  i  =  i  .jo 
8  Pl< I  .  1 )  =  PL ( 3,  1  ) 

KP  =  2 

jo  =  moon 

RETURN 
5  K  =  K  +  T 
XPIK )  =  X 
YP(K)  =  Y 
GO  TO  4 


SUBROUTINE  GRID 

C  THIS  ROUTINE  PREPARES  THE  GRAPH  LIMITS  FOR  PLOTTING  ISOINTENSITY  LINES 
DIMENSION  TITLfc(8),XT(9>,YT{9),PX<3),PY<100),PL<3,100>,FC<20) 

1 , IP2I20) 

COMMON/MCG/  .10,  I  PL  »PX,PL  ,PY,FC,XL  ,  XR , YR , Y? , I P ♦  I  P2 

DATA  TITLF/3-1H  ,  6HPL  UMF  ,  6HC0NT  00 , 2  HRS  /  ,  X  T  ,'4*1  H  ,6HX  (CM)/, 

1  YT/A=?1H  ,6HY  (CM)/ 

CAM.  AlCRT3(n,0,X,Y,2,3,l,l,A2,TITLE,XT,YT,l,l,16.,16.,2,XL,XR,2, 

1 YR ,Y2  ) 

I  PL  =  1 

RETURN 

END 
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SURRUUTINE  C  F  (  X,Y,F,C  ) 

C  THIS  ROUT  I NE  PERFORMS  A  LEAST-SOUARE  CURVE  FIT  TO  9  (X,Y>  POINTS  OF 
C  FUNCTION  F,  WHFRE 

C  F(FIT)  =  cmx=»+?  +  C  ( 2  )  X  +  C(3)Y$*2  C(A)Y  +  C  ( t> )  XY  +  C(6) 

DIMENSION  X(3),Y!3),F(3,3),C(ft),A(ft,6),RT(18),IT<18) 


DO  1  l 

= 

1  ,ft 

Cll)  = 

0. 

DO  1  ,1 

= 

1  ,ft 

1 

A  (  I  ,  J  ) 

n. 

no  ?  i 

1  ,3 

X?  =  XII) 

$s=2 

All  ,1  ) 

= 

A ( 1  ,1)  +  X2*X2 

All  ,2  ) 

= 

A ( 1 ,2  )  +  X2*X(  I  ) 

A  (  ?  ,  ft  ) 

= 

A ( 2  ,ft)  +  XI  1  ) 

2 

A  (  2 , 7  ) 

= 

A  (2 ,2  )  +  X2 

DO  3  I 

= 

1,3 

Y?  =  Y(I) 

*=*2 

A  (  3 , 3  ) 

= 

A ( 3 , 3 )  +  Y2*Y2 

A  (  3  ,  A  ) 

= 

A ( 3 , * )  +  Y2*Y( I ) 

A  (  A  ,  6  ) 

= 

A (A  ,6)  +  Ylll 

3 

A(A,A ) 

- 

MA,h)  +  Y2 

All  ,1  ) 

= 

3.*A<1  ,1  ) 

All  ,2  ) 

= 

3. -A! 1,2) 

A  (  2 , 2  ) 

- 

3.AAI?  ,?  ) 

A  (  3 , 3  ) 

3.*A( 3,3) 

A(  3, A) 

3 • *A ( 3  » A ) 

A  (  A  ,  A  ) 

= 

3.  *A  (  A,  A  ) 

A<  A, A) 

- 

3  .  *  A  ( A  ,  ft  ) 

A  (2 ,6  ) 

= 

3.*A(2,ft) 

A  (  2 , 1  ) 

= 

All  ,2) 

A  (  A ,  3  ) 

= 

A  (  3  ,  A  ) 

DO  A  I 

= 

1  ,3 

DO  A  J 

= 

1,3 

X2  =  XI  I  )**? 

Y2  =  Y ( J ) **2 
XV  s  X(  I  )*Y( J  I 
All, 3)  =  All ,3)  +  X2*Y2 
All, A)  =  All, A)  +  X2*Y(J) 
All, 5)  =  All, 5)  «■  X?  *XY 
A ( ? , 3 1  =  A ( 2  , 3 )  i  XY*Y ( J ) 
A ( 2 , A  I  =  A ( 2 , A  )  +  XY 
A  (  3 , 5  )  =  A ( 3 , 5  >  +  XY*Y2 
Cl’  )  =  Cll)  +  FI  I  ,,1)3X2 
Cl?)  =  C ! 2 )  +  F ( I , J ) *X (  1  ) 
Cl  3)  =  Cl  3)  +  F(1,J>*Y2 
CIA)  =  CIA)  +  H  I ,  J  )  x-Y  (  J  ) 
C  l  ft  )  =  C  (  6  )  +  F  ( I  ,  I ) 

A  Cl 5)  =  Cl  5)  +  F(I,J)*XY 
A  <  3, l  )  =  All, 3) 

A l 3 ,2  )  =  A ( 2 ,3 ) 

A l A , 1 )  =  All ,A) 

A  l  A  ,  2  )  =  A  ( 2  ,  A  ) 

A  ( A ,  5  )  =  A  ( 2 , 3  ) 

A ( 5, A)  =  A ( 2 , 3 ) 

A (5,1  )  =  All, 5) 

A  <  5 , 2  )  =  All, A) 

A ( 2 , 5  )  =  All, A) 

A ( 5 , 3 )  =  A l 3 ,5) 

A ( A,6  )  =  9. 


FUNCTION  ERF (FT  A) 

SEE  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS  WITH  FORMULAS,  GRAPHS  AMO 
MATHEMATICAL  TABLES.  MRS ,  APPLIED  MATHEMATICAL  SERIES  55, 

PAGF  2  99 ,  SECTION  7.1 .2  6  (1964) 

DATA  P/  0.3275911/ 

DATA  A 1 , A? , A  3 , A4 , A  5  /  0.2  5482  9592,-0.2  84496736,1  .421413741, 

1-1  ..453152027,1  .06140543  / 

T  »  1.0/11.  +  P-ARS (ETA ) ) 

S  "  T  =>  ( A 1  +  T*(A2  +  T*(A3  +  T  *  (  A  4  +  A5*T)))) 

ERF  =  1.  -  S:?F XP (  -ET  A *$2  ) 

IF  (ETA.LT.  0.0)  ERF  =  -FRF 

RFTURN 

END 
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FUNC.T  I  ON  TAPU(XX,X,FX,Mn,L,JFX)  .  0010 

D  I  MENS  !  ON  V(MO),FX(MO)  0027' 

M=MO  0030 

If-  IXX-VIHI  1  ,1,3  0040 

1  IF  ( XX-XI ] ) )  2,4,6  003' 

2  T  ARI. E  =  F  X (  1  )  0060 

GO  TO  4  00 7 ( 

3  TARLF=FX(M)  OOKT 

4  WK  I  T  F(  6, 3  )  009( 

6  FORMAT!  iM0,l  OX,  ERROR,  RANGE  OF  TABLE  EXCEEDED  )  0100 

WRITE  (6,101 )  XX,X(1 ),X(M).MO 
101  FORMAT  (  3  E 1 3.6,  110) 

IF  ( IEX.FO.0)  GO  TO  100  0116 

STOP  0126 

6  IF  (  IPX.!  1  0)  GO  TO  14 
Ml  =  1 

M3  =  M  0147, 

M 2  =  M/2  0150 

7  IF  ( X(M2 )-XX)  9,8,8  0166 

P  m 3=  ft?  0170 

GO  TO  10  0180 

9  Ml = M2  0196 

10  J--(Ml  +  M3  )/?.  0200 

IF  (J.EO.M? )  GO  TU  11  021 0 

M2=J  0220 

GO  TO  7  -  0230 

11  IF  (X(J)-XX)  12,12,18  0240 

12  M3=J+1  0?5° 

Ml n  J  0260 

GO  TO  14  0270 

13  M3=J  028C 

M1=J-1  0290 

14  IF  il.NE.O)  GO  TO  15  030' 

TARIE  =  FX(MI  )  +  ( XX-XIM1  ) ) *( FX ( M3  >-FX ( Ml  ) > / < X { M3 )-X  (Ml  )  )  0310 

GO  TO  100  0320 

15  IF  (Ml.GT.l)  GO  TO  16  0330 

11=1  0340 

12=2  0350 

13=3  036f 

GO  TO  17  0370 

16  IF  (M.GT.M3)  GO  TO  18  0380, 

I3=M  0390 

12  =  M-1  0400 

I  1 =M-2  0410 

17  0F32=FX( 13 )-FX( 12 )  0420 

()F21  =  FX<  12  )-FX{  m  043c 

OENOM=(  X(T3)-x(I2n=X(X(13)-X(Il>)*(XU2)-X(Il))  0440 

A  =  (  (  X<  I?)-X(  I!  ))*0F32-(  X(  I  3 » -X  (  12)  )*0F21  l/DENOM  0450 

R=  (  (  X(  12  )##?-/(  1 1  )**2 )*DF32-(X(  I3)#*2-X(  1 2 ) #*2  > *DF2 1 )  /DENOM  0460 

TABLF  =  FX(  I 3)-A*( X( 13 )**2-XX*4?  )+B*<  X ( I  3 )-XX )  0470 

GO  TO  100  0430 

1 H  I 4=M3+1  0490 

1 3=M3  050(' 

I2=M1  0510 

I 1 =M 1 — 1  0520 

DF43=FX( T4)-FX(J3)  0530 

OF32=FX( I3)-FX( 12)  0540 

i)F2 1  =  F X ( 12  1-FX (  U)  0550 
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(A]  have  only  s  <  8  significant  digits,  the  user  may  compute  sigrev, 
the  corrected  number  of  significant  digits  in  the  inverse  by  the 
formula 

sigrev  =  sigdig  -  8.  3  +  s 

To  get  some  idea  how  accurate  sigrev  itself  was,  the  difference 
between  sigrev  and  sigact,  the  actual  average  number  of  significant 
digits,  was  computed  for  many  matrices  of  random  numbers. 

When  alsmo  all  the  elements  of  the  matrix  were  of  the  same  order 
of  magnitude,  the  results  (for  a  batch  of  120  matrices)  were: 

-1.  17  <  (  sigrev  -  sigact  )  <  0.  75 

However,  when  the  elements  of  the  matrix  ranged  over  three  orders 
of  magnitude,  the  results  for  a  batch  of  160  matrices  were: 

-1.  77  <  (  sigrev  -  sigact  )  <  0.  78 

Both  batches  contained  100  5  by  5  matrices,  the  other  matrices 
all  being  10  by  10. 

A  DOUBLE  PRECISION  substitute  for  this  subprogram  is  available. 
I.  Error  Indications  -  ierror  is  an  INTEGER  variable  which  is  used  to 
flag  errors. 

ierror  =  -1  indicates  both  of  the  following: 

1.  Neither  [Al“ ^  nor  [X)  could  be  computed. 

2.  (A)  was  either  singular*  or  very  nearly  so  (usually  the 
former. 

ierror  =  1  indicates  that  no  errors  were  detected. 

Due  to  round-off  errors,  a  singular  matrix  will  rarely  be 
flagged  as  singular  (by  ierror  =  -1)  unless  it  has  a  row  or 
column  of  zeros.  However,  the  value  of  sigdig  will  seldom  be 
greater  than  1.  5  for  singular  matrices. 

Division  by  zero  cannot  occur  in  this  subprogram. 

No  test  for  floating-point  overflow  is  made  in  this  subprogram. 

*When  [Al  is  singular,  a  system  of  simultaneous  equations  [AljXj  -  [B] 
will  have  a  solution  if  and  only  if  [Bl  can  be  expressed  as  a  linear 
combination  of  the  columns  of  [A], 
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Appendix  C 

AUXILIARY  LIBRARY  ROUTINES 

Operating  descriptions  and  listings  are  presented  here  for  LSQP0L  and  SID-- 
two  subroutines  from  the  MDAC  Library.  These  were  not  programmed  under 
this  contract,  but  are  used  by  the  FLAME  code. 

C.  1  LSQPQ5L 

Least  Squares  Polynomial  SUBROUTINE 

A.  Description  -  This  subroutine  subprogram  computes  the  coefficients 
of  the  polynomial  which  best  fits  (in  the  least  squares  sense)  a 
given  function  of  one  independent  variable. 

B.  Use  -  CALL  LSQPOL  (x,  y,  w,  r,  n,  s,  m,  c)  where: 

x  is  the  REAL  array,  dimensioned  at  least  n,  of  the  values 
of  the  independent  variable  of  the  given  function, 
y  in  the  REAL  array,  dimensioned  at  least  n,  of  the  values 
of  the  dependent  variable  of  the  given  function, 
w  is  the  REAL  array,  dimensioned  at  least  n,  of  the  least 
squares  weighting  for  each  point  of  the  given  function. 

(If  an  unweighted  least  squares  fit  is  desired,  each  ele¬ 
ment  of  this  array  must  be  a  1.  ) 
r  is  the  REAL  array,  dimensioned  at  least  n,  which  this 
subprogram  sets  equal  to  the  residuals  at  each  point  of 
the  given  function. 

n  is  an  INTEGER  variable  (or  constant)  denoting  the  number 
of  points  defining  the  given  function.  This  number  must 
be  <  50. 

s  is  the  REAL  variable  which  this  subprogram  sets  equal  to 
the  sum  of  the  squares  of  the  residuals, 
m  is  the  INTEGER  variable  (or  constant)  denoting  how  many 
coefficients  the  least  square  polynomial  will  have.  That 
is,  it  is  the  degree  plus  one  of  the  least  squares  poly¬ 
nomial  desired.  This. number  must  be  input  <  10. 


c  is  the  REAL  array,  dimensioned  at  le?-.>t  nr>  who,  e 
elements  are  set  cqu.il  to  t  .c-  coefficients  of  the  least 
squares  polynomial  beginning  w.Ln  the  conscait  term. 

C.  Support  Level  -  Supported  by  Programming  Systems  and  Support 
Branch  of  Information  Processing  Systems. 

D.  Language  Used  -  FORTRAN. 

E.  Availability  -  On  FOP.TRAN  library. 

F.  Extent  -  4052g  words. 

G  Timing  -  Not  available. 

H.  Restrictions  -  None. 

I.  Ev rot  Indications  -  The  least  squares  problem  is  mathematically 
inconsistent  .tless  m  <  n,  therefore,  if  m  >  n  the  following  note 
is  printed  and  execution  of  the  program  is  terminated. 

LSQPOL  ERROR  -  MORE  COEFFICIENTS  THAN  POINTS 

J.  Method  -  The  method  used  is  completely  described  in  "Polynomial 
Curve  Fitting  With  Constraint"  by  J.  E.  L.  Peck  in  SIAM  Review, 
volume  4,  number  2,  April,  1962. 


This  subprogram  finds  the  Cj's  such  that: 


n 


=  I  wkr 


2  . 

k  is  a  minimum. 


krl 


m 

v  i-1 

rk  =  *k  -  z  ci  *k 

i=  1 

K.  Examples  -  None. 


C.  2  SID 

Single  Precision  Simultaneous  Equation  Solution  and  Matrix  Inversion 
SUBROUTINE 

A.  Description  -  This  subroutine  subprogram  will  invert  a  REAL 

nonsingular  square  matrix,  [A),  and  if  desired,  also  solve  linear 
system(s)  of  simultaneous  equations,  [A]  [X]  =  [B],  where  [B| 
may  have  any  number  of  columns.  Every  column  of  [B]  must  have 
at  least  one;  non-zero  element.  The  dimensions  of  [A]  and  fB] 
are  limited  only  by  the  available  core  storage. 
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Use  -  CALL  SID  (amat,  nrow,  maxrow,  mxcola,  bmat,  ncolb, 
mxcolb,  sigdig,  ierror,  rtemp,  itemp,  scaleb)  where: 

1 .  Input 

amat  is  a  REAL  two  dimensional  array  which  con¬ 
tains  [A].  [A]  will  be  replaced  by  (Al“l  during 

the  execution  of  this  subprogram. 

nrow  is  an  INTEGER  variable  or  constant  which 

denotes  the  number  of  rows  in  [A],  nrow  may 
be  as  small  as  one  (1)  or  as  large  as  available 
core  storage  will  permit. 

maxrow  is  an  INTEGER  variable  constant  which 

denotes  the  maximum  number  of  rows  which 
may  be  stored  in  the  amat  array.  Note  that  a 
matrix  may  have  fewer  rows  and/or  columns 
than  the  array  which  contains  it. 

mxcola  is  an  INTEGER  variable  or  constant  which 

denotes  the  maximum  number  of  columns  which 
may  be  stored  in  the  amat  array. 

bmat  is  a  Real  array  which  contains  [B]  if  [B]  is 

present.  If  [Bl  is  present,  the  bmat  array  must 
have  exactly  maxrow  rows.  [B]  will  be  replaced 
by  [X]  during  the  execution  of  this  subprogram. 

If  no  [B]  is  present,  bmat  may  be  any  variable 
or  constant  of  any  type. 

ncolb  is  an  INTEGER  variable  or  constant  which 

denotes  the  number  of  columns  in  [B] .  If  there 
is  no  [B],  use  ncolb  =  0.  ncolb  may  be  as 
small  as  zero  or  as  large  as  available  core 
storage  will  permit. 

mxcolb  is  an  INTEGER  variable  or  constant  which 

denotes  the  maximum  number  of  columns  which 
may  be  stored  in  the  bmat  array.  If  no  |B]  is 
present,  mxcolb  may  have  any  value. 

2.  Output 

amat  will  contain  lAl"*. 

bmat  will  contain  (X)  if  [B]  was  present. 


sigdig  is  a  REAL  variable  which  will  be  set  equal  to 
an  estimate  of  the  number  of  significant  digits 
in  the  elements  of  the  inverse.  However,  this 
estimate  is  based  on  the  assumption  that  the 
elements  of  [A]  are  all  accurate  to  eight  signi¬ 
ficant  digits.  The  Restrictions  section  describes 
a  simple  adjustment  of  sigdig  which  the  user 
must  make  when  the  elements  of  [A]  have  fewer 
than  eight  significant  digits. 

ierror  is  an  INTEGER  variable  which  is  used  to  flag 
errors.  Its  meaning  is  explained  in  Error 
Indications. 

3.  Intermediate  -  The  following  three  arrays  are  used  during 

calculation  for  the  temporary  storage  of  intermediate  results. 

It  makes  no  difference  whether  the  arrays  are  one,  two,  or 
three  dimensional.  The  three  arrays  need  not  all  be  of  the 
same  dimension. 

rtemp  is  a  REAL  array  which  contains  at  least 
3*  nrow  elements. 

itemp  is  an  INTEGER  array  which  has  at  least 
3’ nrow  elements. 

scaleb  is  a  REAL  array  which  contains  at  least 
ncolb  elements. 

C.  Support  Level  -  Supported  by  Programming  Systems  and  Support 
Branch  of  Information  Processing  Systems. 

D.  Language  Used  -  MSSD  Standard  FORTRAN. 

E.  Availability  -  On  the  FORTRAN  library. 

F.  Extent  -  766  words,  19  words  for  ALOGIO,  and  46  words  for  ALOG. 

G.  Timing  -  Not  available. 

H.  Restrictions  -  This  subprogram  cannot  find  a  solution  to  the  system 
of  simultaneous  equations  [A]  [X]  =  [0], 

An  estimate  of  the  average  number  of  significant  digits  in  the 
elements  of  the  inverse  is  given  by  the  argument  sigdig.  (See  the 
Use  section.)  However,  this  estimate  is  only  valid  when  the  ele¬ 
ments  of  A  have  eight  significant  digits.  When  the  elements  of 
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[A]  have  only  s  <  8  significant  digits,  the  user  may  compute  sigrev, 
the  corrected  number  of  significant  digits  in  the  inverse  by  the 
formula 

sigrev  =  sigdig  -  8.  3  +  s 

To  get  some  idea  how  accurate  sigrev  itself  was,  the  difference 
between  sigrev  and  sigact,  the  actual  average  number  of  significant 
digits,  was  computed  for  many  matrices  of  random  numbers. 

When  alsmo  all  the  elements  of  the  matrix  were  of  the  same  order 
of  magnitude,  the  results  (for  a  batch  of  120  matrices)  were: 

-1.  17  <  (  sigrev  -  sigact  )  <  0.  75 

However,  when  the  elements  of  the  matrix  ranged  over  three  orders 
of  magnitude,  the  results  for  a  batch  of  160  matrices  were:. 

-1.  77  <  (  sigrev  -  sigact  )  <  0.  78 

Both  batches  contained  100  5  by  5  matrices,  the  other  matrices 
all  being  10  by  10. 

A  DOUBLE  PRECISION  substitute  for  this  subprogram  is  available. 
I.  Error  Indications  -  ierror  is  an  INTEGER  variable  which  is  used  to 
flag  errors. 

ierror  =  -1  indicates  both  of  the  following: 

1.  Neither  [ A] ”  ^  nor  [X]  could  be  computed. 

2.  (A]  was  either  singular*  or  very  nearly  so  (visually  the 
former. 

ierror  =  1  indicates  that  no  errors  were  detected. 

Due  to  round-off  errors,  a  singular  matrix  will  rarely  be 
flagged  as  singular  (by  ierror  =  -1)  unless  it  has  a  row  or 
column  of  zeros.  However,  the  value  of  sigdig  will  seldom  be 
greater  than  1.  5  for  singular  matrices. 

Division  by  zero  cannot  occur  in  this  subprogram. 

No  test  for  floating-point  overflow  is  made  in  this  subprogram. 

*When  [A]  is  singular,  a  system  of  simultaneous  equations  (AllXl  -  fB] 
will  have  a  solution  if  and  only  if  (B)  can  be  expressed  as  a  linear 
combination  of  the  columns  of  f A]. 
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J.  Method  -  The  numerical  method  used  is  basically  the  Gauss- Jordan 
elimination  with  selection  of  maximum  pivotal  elements  (full 
pivoting).  A  description  of  the  Gauss-Jordan  method  is  contained 
in  Numerical  Analysis  by  Kaiser  S.  Kunz;  McGraw-Hill,  1957. 

Special  procedures  are  present  to  improve  accuracy  and  also  to 
minimize  the  frequency  of  overflow  and  underflow  when  the  magni¬ 
tudes  of  any  of  the  elements  in  the  A  or  B  matrices  are  either 
large  or  small. 

K.  Example  -  The  coding  below  will  invert  a  matrix,  solve  a  system 
of  simultaneous  equations,  and  print  a  warning  note  if  the  inverse 
has  less  than  about  one  significant  digit.  The  elements  of  the 
matrix  originally  in  AMAT  are  assumed  to  have  an  average  accuracy 
of  eight  significant  digits. 


DIMENSION  AMAT  (10,  10),  BMAT(  10,  1  ),  RTEMP(10,3), 
ITEMP(  10,  3),  1  SCALEB(l) 


NR0W  =  6 


CALL  SID  (AMAT,  NROW,  10,  10,  BMAT,  1,  1,  SIGDIG,  IERROR, 
1  RTEMP,  ITEMP,  SCALEB) 

IF  (SIGDIG  .LT.  2.  )  WRITE  (6,  180)  SIGDIG 
180  FORMAT  (40H 1  WARNING  -  MATRIX  INVERSE  HAS  ONLY  ABOUT, 

1  F5.  1,  2 OH  SIGNIFICANT  DIGITS.  ) 
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c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


SUBROUTINE  SID  (A,  N,  NDRflW »  NDCOLA,  B,  M,  NDCOL°.,  SIGDIG,  I  ERROR , 
.  PIVOT,  INDEX,  SCALER  ) 

SID  -  A  SINGLE  PRECISION  SIMULTANEOUS  EQUATION  SOLVER,  INVERSE 
FINDER,  AND  DETERMINANT  SUBROUTINE 

DIMENSION  AINDROW, NDCOLA!  ,  B I NDRDW , NDCOLR ) ,  PIVOT(N,3), 

.  SCALFB(M),  INDEX (N, 3) 

DOUBLE  PRECISION  DBIGP2 
DATA  DBIGP2 

.  /  737B6976294B3829.04 

EPS  =  l.F-3 
712  EPS  =  EPS/2. 

EPS  P 1 5=  EPS  +  1.5 

IF  (EPSP15  .NF.  1.5)  GO  TO  712 

SIGMCH  =  ALOGlOIl. 522/EPS) 

RIGPW2  =  DRIGP2 
PI  VOTl 1 ,1 )  =  0. 

SCALE  ROWS 

DO  3H  1=1, N 
ROWMX  =  0. 

DO  2  B  J  =  1  ,  N 

IF  (  (  ABS( A( I,J  ))  )  .GT.  ROWMX)  ROWMX  =  ABS(A(1,J>) 

2  B  CONTINUE 

IE  (  ROWMX)  29,  750,  29 
29  CONTINUE 

RONMXI  =  1.  /  ROWMX 
i)0  3?  J  =  1,N 
A  I J  =  At  I ,J ) 

At  I  ,  J  )  =  t  A (  I  ,  J )  *  ROWMX I  )  *  BIGPW2 

IF  ( A (  I  ,  J  )  .EO.  0.)  At  I , J )  =  ( A  I J  *  BIGPW2)  *  ROWMXI 
32  CONTINUE 

IF  (M)  34,  3H,  34 

34  DO  36  J  =  1  ,M 
BIJ  =  R  (  I  ,  J  ) 

B < I , J )  =  ( B ( 1 , J )  *  ROWMXI  )  *  BIGPW2 

IF  ( R ( I  ,  J  )  .EQ.  0.)  B (  I  ,  J  )  =  (BIJ  *  BIGPW2)  *  ROWMXI 
36  CONT INUF 
3B  PIVOTt 1,2 )  =  ROWMXI 

SCALE  COLUMNS 

DO  10  J  =  1  ,  N 
COLMX  =  0. 

DO  4  1=1, N 

IF  <ARS(A< I,J)).GT.  COLMX)  COLMX  =  ABS(A(I,J1) 

4  CONTINUE 

IF  (  COLMX  )  5,  750,  5 

5  CONTINUE 

COLMX!  =  1. /COLMX 
DO  8  1  =  1,  N 
A I J  =  A  (  I  ,  J  ) 

A  ( 1  ,  J  )  =  ( A ( I , J )  *  COLMX I ) *  RIGPW2 

IF  (At  I , J )  .60.  0.)  A ( I , J )  =  ( A I J  *  BIGPW2 )  *  COLMXI 
8  CONTINUE 

10  PIVOT ( J , 3 )  =  RIGPW2  *  COLMXI 


m 
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~  <si-- 


"Tri  !. 


if§!l 


ooo  ooo  ooo 


IF  ( M )  14,24,14 
i  a  oo  ??  .1  =  1  ,M 
COLMX  =  o. 

00  16  I =  1  ,  N 

IF  (ARSIR! I,J) l.GT.  COLMX  )  COLMX  =  ABS<RII,J)> 

1 6  C0NT I NUE 

IF  (COLMX  )  17,  22,  17 

1 7  CUNT  1 00F 

SCALER(J)  =  COLMX  /  RIGPW2 
C0LMX1  =  1  ./COLMX 
DO  20  I--1  ,M 
R1J  -  R  (  I  ,.l  ) 

R!I,J)  =  (  B  (  1  ,  J  )  ♦  COLMX!)  *  R1GPW2 

IF  (  B  {  1  ,  .1 1  .Ft).  0.)  B  (  1  ,  J  )  =  (  R  I  J  -  BIGPW2)  *  COLMXI 
20  CONTINUE 
22  CONTINUF 
2 A  CONTI NUF 

INITlALIZATION 

°MONF  =  1  . 

DO  A2  J=1 ,N 
PIVOT  (.1,1  )  -  0. 

A?  I NDE  X( J , 3 )  =  0 
DO  550  1  =  1  ,N 

SEARCH  FOR  PIVOT  ELEMENT 

ARPIVI=0. 

A 5  DO  1  05  .1=1  ,N 

50  IF  (  INDEX!.),  3  ) -1  )  60,105,60 
60  DO  1 00  K  =1  ,N 

70  IF  ( I NDE  X ( K , 3 ) -1  )  80,100, BO 

80  IF  (ARS(A(J,K)1  -ARPIVI)  100,100,85 

85  I  ROW=.l 

<J0  ICULHM=K 

ARPIVI=ABS(A(J,K  )  ) 

100  CONTINUE 
105  CONTINUE 

IF  (1-1)  115,120,115 
115  IF  (  ARPIVI  .GE,  PI VM IN  )  GO  TO  123 
120  P I VM  IN  =  AR P I V  I 

IF  (ARPIVI )  123,750,123 
123  CONTINUE 

INDEX! ICOLUM, 3 1=1 
PI  VCITI  =A  (  mow, ICOLUM) 

PIVOT! I ,1 )  =  PIVOT  I 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

130  IF  ( IROW-ICOLUM)  140,  260,  1A0 
1  AO  PMONF=— PMONE 
1 50  DO  200  L=1 »N 
160  SWA  P=A 1  I  ROW ,L ) 

170  A( IROW,L)=A( ICOLUM, L ) 

200  A!  ICOLUM, L  )=SWAP 
205  IF  (M)  260,  260,  210 
210  DO  2  50  L  =1  ,  M 
220  SWAP  =  8 ( IRQW»L ) 
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OOO  UOO  o  o  <->  o  o  o  o  o 


2 3n  niiRnw.u  =  r(icolum,l) 

?W  fll  ICHUIM.L  )  =  SWAP 
260  I  NOG  > (  I ,1  )  =  I  ROW 
?  70  INDEX! 1 ,2 )  =  l COLON 

DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

PIVINV=1.0/PIV0T1 
330  A  (  1 C OLUM  ,  ICOLtIM  )  =  R1GPW2 
3 AO  00  360  L=1 .N 

350  A(  ICOUM.L  )  =  A  (  ICOLtIM  ,L  )  *P  1  V  I  NV 
355  IE  (MI  330,  3H0,  360 
360  DO  3  70  L  =1 , M 

370  R ( I COLUM, L )  =  R ( I COLUM, L ) *P I  V  I NV 

REDUCE  NON-PIVOT  ROWS 

380  DO  550  LI =1 ,N 

390  IFILI-ICOLUM)  AOO,  550,  AOO 

400  T=A(L1 , ICOLUM) 

IF  (T)  420,550,420 
420  A'Ll, I C  OL  UM ) =0 .0 
430  DO  450  L=1,N 

450  A(‘.1,L)=A(L1,L)-A(IC0LUM,L)*T 
455  1E(M)  550,  550,  460 
460  DO  500  L=1 ,M 

500  R  (  L  i  »L  )  =  8  (  L  1  ,  L  )  -  R<ICnLUM,L)*T 
550  CONTINUE 

interchange  columns 

600  DO  710  1=1 ,N 
610  L  =  N+ 1 -  I 

620  IF  ( IN0SX(L,1 )-IN0EX(L,2 ) )  630,  710,  630 
630  JROW  =  I  Nf)E  X  (  L  ,  1  ) 

640  JCOLUM= 1  NOE X ( L ,2 ) 

650  DO  705  K  =1  ,N 
660  SWAP=A(K, JROW) 

670  AIK, JROW ) =A ( K , JCOLUM ) 

700  A(K,.!COLUM)  =  SWAP 
705  CONTINUE 
710  CONTINUE 

PIVOT (1,1)  =  PIVOT! 1,1)  *  PMONE 

SIGDIG  =  SIGMCH  -  AL0G10I8IGPW2/PIVMIN) 
IF  (SIGDIG  ,LT ,  .85)  SIGDIG  =  0. 

UNSCALE  INVERSE  AND  SOLUTIONIS) 

DO  720  J=1 ,N 
ROWMXI  =  P I VOT ( J , 2 ) 

00  720  1=1, N 

IF  tROWMXI  .LT.  1.)  GO  TO  715 
A (  I  ,  J  )  =  ( A ( I , J  )  *  PIVOT (1,3))  #  ROWMXI 
GO  TO  720 

715  A ( I  ,  J  )  =  A (  I  ,  J  )  *  ( PIVOT (1,3)  *  ROWMXI) 
720  CONTINUE 

IF  (M)  725,  735,  725 


725  1)0  730  .1  =  1 

ROHM  X I  =  SCALBUJ) 
no  730  1=1. N 

.«™«i 

7?fi  »’u"™*n.0l  •  .PIVOT.. ,31  *  HOWMX I 

730  CONTINUE 
735  CONTINUE 

I  ERROR  =  1 
RETURN 

750  I  ERROR  =  “1 
SI  GO  10  =  0. 

RETURN 


SUARPOTINF  I.SOPOL  (  X  1  ,  Y 1  ,  W!  ,  RD  ,  N  ,  SUM  ,M  ,COEF  ) 

11  IMF  NS  ION  XI  ( 1  )  ,Y1  ( l  )  tWI  ( 1  )  ,RO(  1  )  »COFF(  1  ) 

COMMON  /  1/ I.SOP/ 

.  X(50)  ,Y(50) ,W( 50) tB(ll l»AC 11 ), S(ll >  tDSQIll ) 

.  ,P(50),P0(50),C<11 > 

OOURLF  PRECISION  X  ,  Y ,  W ,  B  ,  A  ,S  ♦  PO  ,  P  ,DSQ,C  .  G 
DO  1 00  I  -  1  i  N 

xu  )  =  xi  <  i ) 

Y (  I  )  =  YI (I) 

100  W (  I  )  =  Will) 

CALL  ORPOLG ( M-l *  X  *  Y, W, N» B  »  A , S , DSQ  t  PO  « P , G, E ) 

I  F ( F  )  1 1 00,1 50, 11 00 
150  CALL  ORPOLC I M-l ,C , S , A ,R ,G, P, PO ) 

DO  2  00  1  =  1 ,M 

200  COfcFI  1)  =  cm 

SOM  =  0.0 
DO  500  I  =  1,N 
K  =  M-l 
Z  =  COEF(M) 

IF  (K)  300  » 350 »  300 
300  Z  =  Z*XI ( I  )  +  CGEFtK  1 
K  =  K-l 

IF  (K)  350,350,300 
350  ROIII  =  7  -  YKII 
500  SOM  =  SOM  +  ROU) 

ret  urn 

1100  WRITF  (6,1200) 

1200  FORMAT(43H1LSOPOL  ERROR-MORE  COEFFICIENTS  THAN  POINTS) 
STOP 
END 
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SURROUTI  NE  ORPOLGIN,  X, Y,W,M,  R,A,S,DSQ,P0,  P,  GAMMA,  ERROR 

DOURI.F  PRECISION  X,Y,W,R,A,PO,P  ,DSQ , WPP , WXPP , WY P , TE MP 
DOURM  PRECISION  WPPO, SOM, GAMMA, 1  ,MEAN,S 

EOUl VAL  FNCF  I T  >  TFMP, SUM ) 

DIMENSION  X(1 I.YIl ) ,WI1) ,R(1) ,A(1 ) ,S(1) ,P0(1) ,P  (1), 
l  nsum 

SUM=0. 

DO  6  1=1 ,M 
SUM=SUM+X( I ) 

ME AM=  SOM/ FLOAT (M) 

T  =  0. 

DO  7  I = 1 , M 

T=  AMA  XI ( T , AH  S ( X ( 1 (-MEAN ) ) 
r  TT  =  DARSI XIII-  MEAN) 

IF  (IT  .GT.  T)  1  =  TT 
GAMMA =  T/2 . 

DO  R  T  =1  ,M 

X  (  I  )  = ( X ( I  l-MEAN) /GAMMA 
EPR0R=0. 

NO=M-N-l 

IF ( N0)102 ,103,103 
E  RROR  =  1 . 

GO  TO  105 
B  =0  • 

0SO=0. 

WPP=0. 

DO  106  J  =  1  ,M 
P  (J)=l. 

P0( J )=0. 

WPP=WPP+Wf J ) 

)F(N 0)106, 106,107 


107 

DS0=DSQ+WU>*Y(J)**2 

4  * 

106 

COOT 1N0E 

_ 

NA  =  N+ 1 

DO  109  1=1, NA 

WXPP=0. 

WYP=0. 

DO  110  J  =  1  ,M 

T  EMP  =VJ  ( J  )*P(,I  ) 

I F ( NA- I >111,111,112 

£ 

112 

WXPP=WXPP+TEMP-X(J )*P< J  ) 

* 

111 

I F (  M  -  1)110,113,113 

113 

WYP  =  WYP  +  TEMP*Y( j  ) 

y 

no 

CONTINUE 

& 

I F ( M  -  I >114,115,115 

& 

(I" 

115 

SI  I  )  =  WYP/ WPP 

& 

1 1 A 

I F ( NA-M  >666,117, 117 

§r 

666 

I F ( 1-1 1102,667,116 

% 

667 

DSO  =  DS0-S!**2*WPP 

Jr 

GO  TO  117 

'1| 

116 

DSO< I )=DSCH  I  — 1 >  — S  C I )**2*WPP 

117 

IFINA-I 1109,109,119 

m 

119 

At  I )=WXPP/WPP 

WPPO=WPP 

WPP=0. 

DO  120  .1  =  1, M 

TEMP=( XI J )-A( I ) )*P( J )-R( I )*PO(J ) 
WPP=HPP+W I J ) *T  EMpM*2 


PO(J)  =  P(J) 

120  P ( J )  =  TEMP 

R( 1+1 )=WPP/WPP0 

109  com  i  Niie 

IP  (N)  1  AO,  155,1  AO 
1  AO  00  150  1  =  1  , N 

A  (  I ) =A  (  I  )*GAMt-'A  +  PE  AN 
150  B ( I  )  =  R ( I ) SGAMMA 
155  DO  1 50  I =1  ,M 
160  X(  I  )  sGAMMA-^X (  1  1+MEAN 
105  RETURN 
END 


SURROllT  INF  OR  POL  C  <  N  ,C  ,  S  ,  A  ,  R  ,  G ,  P  ,  PM  ) 
DOUBLE  PRECISION  C , S , A ,R , G , P , PM, T1 , T2 
DIMENSION  C(1),S(1) ,A(1) ,R(1) ,P(1) ,PM(1) 
M  =  AN  1 

00  300  L  =1  , M 
C (L )=0. 

PM  I L ) =0  . 

300  P(L)=0. 

P=1  . 

C  =  S 

IF  IN)  301,303,301 
301  00  302  1  =  1  ,N 
T2=0. 

J  =L  +  1 

DO  302  M  =  1  ,J 

T1  =  (T?-ML)*P(M)-B(L)*PM(M)  1/G 
T2  =  P  (  M  ) 

PKI  M ) =P ( M ) 

P  (  M  >  =T  l 

302  C(M)=C(M)+T1*S(L+1 > 

303  RETURN 
END 
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